J. Plasma Fusion Res. Vol.92, No.6 (2016)４８２‐４８７

講座

今日からはじめるメニーコアアクセラレータ
Let's Start a Numerical Simulation Today Using the Many-Core Accelerator

１．メニーコアアクセラレータを研究室で使ってみよう
1. Let’s Use Many-Core Accerelator in Your Laboratory
八 柳 祐 一１），田 邊 恵 李１），堀 越 将 司２）
YATSUYANAGI Yuichi１）, TANABE Eri１）and HORIKOSHI Masashi２）
１）

静岡大学教育学部，２）インテル株式会社ソフトウェア＆ソリューション統括部（原稿受付：２
０
１
６年３月１
８日）

大学の研究室レベルで保有可能でありながら，マッシブな（大規模シミュレーションが可能な）計算環境を提供してくれる Intel" Xeon PhiTM コプロセッサ（アクセラレータ）の紹介を行う．アクセラレータは GPU（Graphics Processing Unit）と同様に，PCIe バスに挿す形で提供される並列計算機である．本章では，アクセラ
レータの特徴，及び導入の仕方を解説するとともに，非中性プラズマシミュレーションを例として取り上げ，計算可能なシミュレーション規模を視覚的に示す．
Keywords:
parallel computing, OpenMP, parallelize, vectorize, N-body simulation

１．
１ 初めに

ming Interface］）を 用 い て コ ー ド を 作 成 す る 必 要 が あ

１９９０年代の PCと言えば，
CPU は i386
（25 MHz）
や i486DX2

る．すなわち，インテルや AMD 製の CPU を搭載する PC（66 MHz），Pentium（166 MHz），PentiumII（300 MHz），メ

で動く並列コードは，そのままの形では，GPU などで走ら

イン メ モ リ は 16 MB〜256 MB が 当 た り 前 で あ り［１，
２］，

せることはできない＊２．

１０００粒子程度のちょっとした ! 体シミュレーションを行

このように，ハードウェア環境ごとにプログラミングやおうとすると，計算機センターに設置してあるスーパコン

チューニングを行わなくてはならないのは，研究者にとっピュータを使わざるを得ず，長いジョブ待ち行列にうんざて本質的ではない．そこで，当研究室では，PC 用の並列りしたものである．

コードを再コンパイルのみで超並列実行できる計算機，Intel" Xeon PhiTM コプロセッサ（以後， アクセラレータ

その後，大学の１研究室レベルで保有可能なデスクトッ
プスーパコンピュータの時代が，１９９０年代後半〜２０００年代と呼ぶ）に注目した．

に到来する．重力やクーロン力といった相互作用力の計算

本章では，アクセラレータの特徴，および使い方についのみに特化した専用計算機（GRAPE/MDGRAPE シリーて簡単に紹介するとともに，非中性プラズマでのdiocotronズ）により，１０４〜１０５粒子のダイレクト・シミュレーショ

不安定性を例にとり，どの程度の規模のシミュレーション

ン が 現 実 的 な 時 間 で 可 能 に な っ た［３，４］．さ ら に そ のが可能なのか解説をする．具体的な速度比較は次章で扱う

後，２０１０年代にかけてGPUが一般的になり，PC１台分の大（図１）．

きさの計算機で１TFLOPS＊１の計算速度を実現するのも

１．
２ アクセラレータの特徴，使い方

夢ではなくなった．
しかし，これらの専用計算機や GPU では，それに特化しアクセラレータは，PCIeバスに挿すボードとして提供さたプログラミング言語（又は，API［Application Program-

れ，サイズは一般的な GPU と同じである．
６４ビット命令が＊１

浮動小数点数の乗算などを１秒間に 1 T=1012 回計算する能力．

＊２

近年，GPU と PC 上でシームレスに動くような API 環境も整備されつつあるが，性能を引き出すのが困難なため，あまり普及していない．

Faculty of Education, Shizuoka University, Shizuoka City, SHIZUOKA 422-8529, Japancorresponding author’s e-mail: yatsuyanagi.yuichi@shizuoka.ac.jp

４８２!2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

1. Let's Use Many-Core Accerelator in Your Laboratory

Y. Yatsuyanagi et al.

MPSS（後述）：MPSS 3.4.2
以下では，この環境を前提に解説する．
１．アクセラレータを搭載した PC を購入する
アクセラレータには，メモリ容量と コ ア 数 に よ って，５７コア＋６G バイト，６０コア＋８G バイト，
６１コア
＋16 G バイトの３種類の製品が提供されている．アクセラレータを動かすホスト PC には，Xeon プロセッサが必須である＊３．
２．ホストコンピュータに OS をインストールする
図１

アクセラレータを動作させるためのソフトウェア（デPC に設置されたアクセラレータの様子．

バイスドライバ等）は，Manycore Platform Software
実行可能な Pentium ベースの CPUコアが６１個，コアあたり

Stack（MPSS）と呼ばれる．MPSS でサポートされる

４スレッド実行可能で，合計２４４並列の計算が可能である．

OS（Linux，又は Windows）のバージョンには制限がGPU などの並列度を見慣れると，「２４４」並列というのはあるので，それに適合する OS をインストールする．
３．インテル製コンパイラをインストールする少ない印象を受けるが，GPU のストリーミングプロセッサと異なり，アクセラレータの CPU コアは，PC に搭載されアクセラレータ上で，Xeon プロセッサ用にコンパイル

るCPUコアと同様の能力を有する高機能なものである．各

されたバイナリがそのまま動く訳ではなく，アクセラ

コアでは５１２ビット長のベクトル演算が可能で，倍精度浮レータ用にコンパイルし直す必要がある．gcc を用い動小数点数であれば８要素の同時計算が可能である．まることもできるが，現状ではパフォーマンスに問題がた，理論上は１コアあたり２スレッド以上の並列度でピーあるので，インテル製コンパイラを使用するのが良い．
４．MPSS をインストールする

ク性能を達成することが可能である．
GPU とアクセラレータの最も大きな違いは，アクセラ

MPSS に添付される README ファイルに従って作業

レータ内で独立した OS が稼働していることだろう．すなを進める．RPM パッケージの形式でドライバ群が提わち，ホスト PC からは，LAN で接続された計算サーバの供されるので，問題が発生することはほぼ無いと思っように見える．一般的に GPU では，初期化処理などをホスて良い．
５．MPSS の初期設定

トCPUで計算し，相互作用計算などの重たい計算部分のみを GPU にまかせる「オフロード」方式を採用する．このたMPSS のインストールが完了したら，め，メモリバスに比べて低速な PCIe バス越しに大量の配・初期化をして，アクセ ラ レ ー タ 用 の ま っ さ ら な

列データをやり取りせざるを得ず，いかにこの通信を減らLinux 環境をホスト PC 上に構築

すかがGPU上のプログラミングでの腕の見せ所となる．一・ユーザを追加

方，アクセラレータでは，重たい計算部分のみを「外注」す・ホスト PC とアクセラレータとのファイルの受渡し

るオフロード方式に加えて，アクセラレータに ssh などでのための NFS（Network File System）を設定

リモートログインし，アクセラレータ用にコンパイルされなどを，ホスト PC 側から micctrl コマンドを使用したプログラムをシェルから起動する，といった使い方がでて行う．一般的な Linux PC の管理を行った経験があれ

きる．すなわち，イメージとしては，計算サーバに ssh でロ

ば，１５分程度で作業は終了する＊４．なお，シェル上でグインし，そこでプログラムをコンパイルして走らせるイ

ホスト PC 側からアクセラレータにリモートログイン

メージに近い．この場合，計算対象となるデータが PCIe

する際に，いちいちパスワードを入力するのは面倒な

バス上を行き来するオーバヘッドがなくなるため，実行速ので，公開かぎ認証を用いてパスワード無しでログイ

度面での恩恵が大きい．

ンできるようにしておくと便利である．
６．アクセラレータの起動

アクセラレータを使い始めるための手順を簡単に示すと，次のような流れになる．なお，当研究室で使用してい「アクセラレータの起動」とは，アクセラレータ用のるソフトウェア等の名称，およびバージョンは，以下のとOS をアクセラレータ内で起動することを意味する．

おりである．

ホスト OS にサービスとして MPSS が登録されている
ので，OS 環境に応じて＊５，

OS：CentOS 6.6
コンパイラ：インテル Parallel Studio XE 2015 Com-# service mpss start

poser Edition for C++ Linux

または，

＊３

Core プロセッサでは動作が保証されていない．

＊４

NFS の設定の際にホスト PC 側のファイアウォールの設定変更を忘れないように！

＊５"#"は管理者権限（root）のプロンプト，"$"は一般ユーザのプロンプトを表す．

４８３

Journal of Plasma and Fusion Research Vol.92, No.6 June 2016

れる．

# systemctl enable mpss

軸方向に印加された磁場，および電子と導体壁との間のとして，アクセラレータを起動する．概ね１分程度で誘導電場中に置かれた電子の運動方程式は，電子の質量を起動する．

#"，負号をつけた素電荷を !"とすると，７．プログラムの準備#""##!(#"# %!"'#%

OpenMP（どのように並列化するかをコンパイラに指示する規約．次章参照）でコーディングされたソー（１）

スファイル（今は C 言語を仮定する）を準備し，インと書ける．高周波のサイクロトロン運動を落とすため，左テル製の C コンパイラを

辺の時間微分を落とし #について整理すると，#% "#!
"!))$ icc -mmic simulation.c -o simulation.xと起動し，アクセラレータ用のバイナリを生成する．

" $
% #$
!!

コンパイル環境と実行環境が同一の場合，インテルコ
ンパイラは-xHOST をオプションに指定することによ（２）発揮できるオプションを自動設定する．アクセラレーを得る．ここで，円筒容器の軸方向を &軸にとり，磁場を$とした．さらに，$を用いて !%!!$&方向の単位ベクトル $タの場合は，-xHOSTの代わりに，-mmicを指定し，ア

電場を静電ポテンシャルを用いて "%!&#とすると，り，当該システムのCPUのベクトル演算能力を最大限

クセラレータ用のバイナリを生成する．生成されたバ#%

イナリは，アクセラレータと NFS で接続されたディレ
クトリにコピーする．また，アクセラレータ内に存在

$
$
#&#
!!

（３）

しない共有ライブラリもあるので，ldd コマンドを使が得られる．この式の形と２次元流に対する流れ関数の形

用し，足りない共有ライブラリを調べ，それらのファとの類推から，イルもホスト PC からバイナリが置かれているディレ・サイクロトロン運動を無視した電子の流れは，非圧縮

クトリにコピーする．

#%!）である（&$

８．実行

・静電ポテンシャル #と流れ関数 %には，%%#!
!!なるアクセラレータに ssh でリモートログインし，先ほど比例関係が成り立つ

コピーしたファイルを実行する．Linux のシェル上でことがわかる．（３）式の rot をとると，渦度を $&として，普通にプログラムを起動するのとまったく同じである．

$%&##
$&$

$ ./simulation.x
以上が，アクセラレータの導入，およびシミュレーション

%

&"#
$
$
!!

%

"
$$
$
"
!!!

の実行過程の概要である．OpenMP でコーディングされたコードを持っているのであれば，GPU よりも気軽に始められる点を，改めて強調したい．

（４）となる．ここで，最後の式変形で，Poisson 方程式を使っ

１．
３ シミュレーション例

た．$は電子の数密度である．この式から，

ここでは，非中性プラズマに対応した点渦系シミュレー

・渦度 $&は，電子の数密度 $に比例する

ションの紹介を行う．点渦系は，渦度方程式（６）を粒子の集合体で離散化するモデルの一つであり，計算機シミュことが分かる．これにより，電子の輝度分布を測定するこレーションでもよく使われる技法である．本節では，アクとは，すなわち，渦度を測定することと同等となる．さらセラレータを使用すると，どの程度の N 体シミュレーショ

に，電子の数密度は空間中での生成消滅がないので，連続ンが可能になるのか，視覚的に示す．シミュレーション速の式

度に関する定量的な評価は，次章に掲載する．非中性プラ

&
$
"#$
&$%!
&
%

ズマ実験に関する解説は［５］と重なる部分もあるが，読者の便宜のため，特に点渦系との対応に関する部分のみ，こ（５）を満たす．この式の両辺に "!
"'!!!(をかけると，

こで解説する．
１．
３．
１ 非中性プラズマと点渦シミュレーション

&
$&
"#$
&$&%!
&
%

プラズマの電気的中性条件を破ったプラズマが非中性プ
ラズマであり，最も極端な例の一つが電子のみから構成さ（６）れた純電子プラズマである［６］．電子プラズマを閉じ込めなる渦度方程式が得られ，電子の数密度 $の時間発展方程る装置としては，円筒真空容器の軸方向に磁場，軸方向両式は，２次元非圧縮非粘性の Euler 方程式と同等であるこ

端に負電位を印加する Penning-Malmberg トラップが使わとが示される．
４８４

Lecture Note

1. Let's Use Many-Core Accerelator in Your Laboratory

Y. Yatsuyanagi et al.

以上より，サイクロトロン運動を無視した電子の運動を

この自転周期と同じ時間刻みで，"$#!!まで行う．すなわ

数値的に追跡するためには，（６）式を解けば良いことがわち，全シミュレーション時間中に渦塊は８回程度自転運動

かる．このために，点渦法を用いる．点渦法は，渦度をする．
時間発展結果を図２（内側半径＝!
!
!!），図３（内側半$Dirac のデルタ関数 $&!'で離散化する．!
#，##は，#番目

径＝!
!
!!）に示す．図２では，"$#!頃にドーナツが崩れ%の点渦の位置ベクトルと循環（強度）である．

始め，モード３が励起されたのち，"$'!には熱平衡分布
!'
!!!
&&&
$! ##$&#'

を目指し，中央にピークを持つ一山分布へと遷移し始め（７）#

る．その後，少しずつ，軸対称な分布を目指して緩和を続!
#の時間発展方程式は，以下の Biot-Savart積分によって与ける．一方，図３では，図２と同様に，"$#!頃にドーナえられる．これは，（７）式を（６）式に代入すると得られツが崩れ始め，モード６が励起される．これを渦結晶配位

る．

と呼ぶ．その後，すぐに熱平衡分布をめざすのではなく，
しばらく六つの渦塊が対称的な配位を保ったまま公転を続

*
!!!'
&
#"
)!
"
#
$! ! #$ # $ #
)%
%$%
#
!
(
!!
(
#
$
$#

ける．その後，"$(!での右下の渦塊群のように，ある二つの渦塊のすそ野同士の接触事故がおこると，渦塊同士の

*
!!! '
&
#"
"
" ! #$ # $ #
%$
#
!
(
#!!
$(

融合が起こり，６→５→４→…と渦塊数が減少し続け，一（８）

山分布の熱平衡解を目指す．しかし，"$"'!以降，一旦，##

なお，ここで半径 !の円筒境界の効果を !$$! !
!
(
$"
$(

渦結晶配位に比べて軸対称性が高い配位に落ち着くと，緩に置いた鏡像渦で取り入れた．（８）式の和は，全点渦に関和速度が遅くなり，十分緩和が進む前にシミュレーション

してとる．つまり，すべての点渦に対する )!
)%を求める#"は終了する．

ためには点渦数の２乗に比例する時間計算量を必要とする

線型不安定性解析では，壁半径，ドーナツの外側，及びため，１０００粒子を超えるシミュレーションとなると，速い

内側半径から線型不安定なモードを決定することができ計算機を使いたくなってくる．我々は，この Biot-Savart
積分を高速に計算できる計算機として，Intel! Xeon PhiTM
コプロセッサ（アクセラレータ）を選択した．
１．
３．
２ 結果
ここでは，アクセラレータを使用して，どの程度の規模の点渦シミュレーションが可能か，視覚的に提示する．具体的な計算速度等は，次章で解説する．
非中性純電子プラズマで見られる有名な不安定性に，
diocotron 不安定性がある．Kelvin-Helmholtz 不安定性と言った方が通りが良いかもしれない．磁場に垂直な断面内で電子の初期分布がドーナツ形状＊６をしていると，シア流（せん断流）により，ドーナツは幾つかの塊に分裂する．
ドーナツが分裂する塊の個数（線型不安定なモード）など
は線型安定性解析により細かく調べられており，シミュ
レーションが定性的／定量的に正しい計算を行っているかチェックする対象として最適なので，ここでは diocotron
不安定性に関するシミュレーション結果を掲載する．
点渦シミュレーションのパラメータを表１に示す．ドー
ナツの外側半径と１点渦あたりの循環，点渦の数密度を一定にして，ドーナツの内側半径をパラメータとした．時間
発展は，（８）式で追跡する．長さの単位は，壁半径!!で規格化している．特徴的時間スケールについて，24.7 という値は，半径 !
!
!!の円内にドーナツと同じ数密度で点渦を&一様に配置した渦塊の自転周期から求めた．時間発展は，
表１

シミュレーション設定．

ドーナツの外側半径 !"ドーナツの内側半径 !#鏡像点渦を除いた実点渦数
特徴的時間スケール "!
＊６

!!
&!!
!!
$!!，!!
%!!
６
１
４
１，３
３
８
２
24.7

図２

時間発展結果．ドーナツの内側半径＝0.3R0 の場合．

３次元的にはちくわ状である．ちくわを包丁で切った断面がドーナツ形状になる．

４８５

Journal of Plasma and Fusion Research Vol.92, No.6 June 2016

表２

線型不安定なモード．

内側半径 #"

線型不安定なモード

!!
##!
!!
$#!

3，2
5，4，6，3，2

さらに，物理的に正しい結果を得ていることの確認のため，上記時間発展結果に対して，! 関数の時間発展を追跡した結果を図４に示す．有限個の正方格子で点渦の配位空間を区切り，%番目のセルに含まれる点渦数を&%で表すと，! 関数の値は次式で与えられる．
!!! &%(
)'&%

（９）

%

! 関数の一般的性質として，その値は時間と共に広義単調減少し，熱平衡に達した段階で &!"&'!!となる．また，準安定状態などに系がトラップされても，&!"&'!!となる
ため，系が現在どのような状態にあるのか，系の緩和に関する貴重な情報源となる．図２のシミュレーションに対応
する ! 関数の時間発展は，$!$!前まで，ほぼ一定値を保っているが，その後単調減少する．これは $!%!以降の軸対称性が高い分布において，熱平衡状態へ向けた緩和が進んでいることを表している．一方，図３のシミュレー
ションに対応する !関数の時間発展は，$!#!程度で値が減少するが，すぐに一定値を保つようになる．これは，準安定な渦結晶配位に留まっていることを表す．その後，値が減少に転じているのは，渦結晶配位が崩れ，熱平衡状態への緩和が進んでいることを表す［８］．

図３

１．
４ まとめ

時間発展結果．ドーナツの内側半径＝0.4R0 の場合．

本章では，アクセラレータの特徴，および使い方につい．シミュレーション設定に対応した線型不安定なモーて簡単に解説を行うと共に，どの程度の " 体シミュレー

ドを不安定な順に列挙した結果を，表２に示す．表２から，ションが可能であるか，視覚的に示すことを目的として，

図２での最も線形不安定なモードは３であり，この結果を

非中性プラズマ実験に対応した点渦シミュレーションの結＊７

る

見る限り，定量的に正しい計算を行っている証拠の一つとなる．しかし，図３での最も線形不安定なモードは５になるはずだが，シミュレーションではモード６が現れている．このような現象は，初期の点渦配置で，周方向の点渦数が，たまたま６の倍数になっていたりすると現れることがある［７］．乱数を用いて特定のモードにエネルギーを注入しないように注意をしないと，このような問題が発生することを，教育的な例として指摘したい．
図２，３に掲載したシミュレーションの実行時間は，それぞれ５１８分，１８３分である．本シミュレーションでは，境界条件を鏡像渦により表現しているため，Biot-Savart 積分
で加算を行う粒子数は，それぞれ２倍の１２２８２粒子，６７６４粒子である．長年，専用計算機 MDGRAPE シリーズで点渦シ
ミュレーションを行ってきた経験からすると，N 体シミュ
レーション研究に十分耐えるだけの速度が実現されている

図４

と感じる．
＊７

H 関数の時間発展．

Davidson の教科書［６］に掲載されている線型不安定性解析では，電子分布の内側にも導体壁がある，すなわち，２重連結領域内に電子群を閉じ込めるような設定を許容しているが，今回のシミュレーションでは内側には導体壁が存在しないので，上記三つのパラメータのみに依存する，と書いている．

４８６

Lecture Note

1. Let's Use Many-Core Accerelator in Your Laboratory

Y. Yatsuyanagi et al.

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

数万粒子程度の点渦シミュレーションであれば，アクセラ
レータにより十分実行可能であることがわかる．地方国立大学の１研究室レベルで保有可能なマッシブな計算環境として GPU 以外の選択肢が増えたことは良いことだろう．

静岡大学教育学部准教授．主な研究分野：GPU や Xeon Phi などの小さくても速いコンピュータを駆使したシミュレーション研究．特に，長距離相関に支配された N 体系の自己組織化に興味があ
ります．２
０
１
１年７月号小特集の時のように我が子との写真を掲載したかったのですが，このアイデアは子どもたちによりゴミ箱へダンクシュート…．やむなく，ガ○プラと．

!!!!!!!!!!!!!

いち!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!

ほり

こし

まさ

し

インテル株式会社，シニアアプリケーショ
ンエンジニア．大阪大学大学院工学研究科
博士課程修了後，２０
０
６年から現職．これまで HPC の分野でスーパーコンピューター
向けに ISV や製造業，大学，研究所のアプリケーションの並列化，最適化およびベンチマークに従事してきた．

!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

堀 越 将 司

た

なべ

え

り

田 邊 恵 李兵庫県立大学シミュレーション学研究科前期博士課程基盤研究室所属１年．学部の卒論では，Xeon Phi を用いて非中性プラズマの緩和過程に関するシミュレーションを点渦モデルにて行った．今年，遂に母が結婚した歳に自分が達し，就職活動と同レベルの焦りとプレッシャーを感じている．

!!!!!!!!!!!!!!!!!!!!!!!!!!

やつ やなぎ ゆう

八 柳 祐 一!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!

参考文献!!!!!!!!!!!!!

［１］SE 編 集 部：僕 ら の パ ソ コ ン３
０年 史（翔 泳 社，東 京，
２
０
１
０）
．

［２］坪山博貴：CPU の謎（ソーテック，２００５）．
［３］伊藤智義：スーパーコンピューターを２０万円で 創 る（集英社（新書）
，２
０
０
７）
．
［４］杉本大一郎（編）：専用計算機によるシミュレーション（朝倉書店，１９９
４）．
［５］際本泰士：プラズマ・核融合学会誌 77, 338 (2001) .
［６］R.C. Davidson, Physics of Nonneutral Plasmas (AddisonWesley, Redwood city, CA, 1990).
［７］Y. Yatsuyanagi et al., Phys. Plasmas 10, 3188 (2003).
［８］Y. Yatsuyanagi et al., J. Phys. Soc. Jpn. 84, 014402 (2015).

果を例示した．これから，境界を表現する鏡像粒子込みで!!!!!!!!!!!!!!!!!!!!!!!!!!

４８７

J. Plasma Fusion Res. Vol.92, No.7 (2016)５
６７‐５７２

講座 今日からはじめるメニーコアアクセラレータ

２．メニーコアアクセラレータ向けにプログラムをチューニングしよう
2. Let’s Optimize a Program for Many-Core Accelerator
田 邊 恵 李１），八 柳 祐 一１），堀 越 将 司２）
TANABE Eri１）, YATSUYANAGI Yuichi１）and HORIKOSHI Masashi２）
１）

静岡大学教育学部，２）インテル株式会社ソフトウェア＆ソリューション統括部（原稿受付：２
０
１
６年３月１
８日）

前章に引き続き，大学の研究室レベルで保有可能でありながら，マッシブな（大規模シミュレーションが可能な）計算環境を提供してくれる Intel" Xeon PhiTM コプロセッサ（アクセラレータ）の紹介を行う．本章では，一般的なシミュレーションコードのチューニング技法として，ベクトル化，並列化といった概念や，並列化に際して必要となる OpenMP の概要を紹介すると共に，アクセラレータでの性能がどのように改善されるのか，具体的な事例を引用しながら定量的な情報提供を行う．
Keywords:
OpenMP,optimization,parallelize,vectorize,N-body simulationに変換することを指す＊１．（１）〜（３）式での加算は スカ

２．
１ 初めに
本章では，Intel" Xeon PhiTM コプロセッサで走るシミュ

ラ命令 と呼ばれるのに対して，（７）式での加算は ベクトル命

レーションプログラムの開発の仕方を，前章で述べた点渦

令 または SIMD（Single Instruction Multiple Data stream）

系を例に説明する［１］．以降，Intel" Xeon PhiTM コプロセッ

命令 と呼ばれる．まとめて処理できる要素数 !をベクト

サをアクセラレータと呼ぶ．

ル長と呼び，計算機（プロセッサ）によって異なる．アクセラレータでは，５
１２ビット長のデータに対応した SIMD

２．
２ 高速化技法

命令があるため，倍精度浮動小数点数（６４ビット長）で８

ここでは，シミュレーションプログラムの高速化に使われ個の演算が１命令で処理できる．すなわち，（１）式のようる代表的技法であるベクトル化，並列化［２］，および，今やな演算を８回行うループ処理があったとすると，スカラ命GPU や CPU などを含めて一般的な高速化技法の一つと令ならば８個の加算命令が必要なのに対して，ベクトル命

なった FMA（Fused Multiply and Add）の利用に分けて概令ならば１個の加算命令で済むため，大幅な高速化が望め説する．なお，並列化に関連する話題として，並列処理をる．しかし，コンパイラは，ループ処理を自動でベクトル行なうための宣言子である OpenMP についても触れる［３］．

化してくれる訳ではない．なぜならば，ベクトルの要素間

２．
２．
１ ベクトル化に依存性がないとコンパイラが判断できない時には，ベクベクトル化とは，今まで，一つずつ順番にトル化されないためである．ベクトルの要素間に依存性が

$
!""
!!#
!

（１）

$
"""
"!#
"

（２）

$
#""
#!#
#

（３）

「ある」とは，例えばベクトル !"#"""%$の要素 "
#
!!
"!
#!
が"
!と "
"の関数となっているような場合，"
!と "
"の値が
決定しない限り，"#の値を決定できないことを表す．このような依存性を落とすためには，プログラムの構造や変数…

定義を根本的に見直す必要があるときもあるので，注意が必要である．

という具合に処理していた演算を，
!"#
"
"
"
%!
!$
!!
"!

（４）

""#
#
#
#
%!
!$
!!
"!

（５）

#"#
$
$
$
%!
!$
!!
"!

（６）

２．
２．
２ 並列化
一般に，単一のオペレーティングシステムで管理される
一つのシステム内での並列化と，別々のオペレーティング
システムで管理される多数のシステム間での並列化では，用いなくてはいけない技 術 が 異なり，前 者では OpenMP，というベクトルに対する一つの演算#"!!"

後者では MPI（Message Passing Interface）が代表的な手（７）法である．今回のアクセラレータは前者に該当するので，＊１ （４）式のようなベクトルデータは，プログラミング言語上では通常，配列で表現されることに注意してほしい．

Faculty of Education, Shizuoka University, Shizuoka City, SHIZUOKA 422-8529, Japancorresponding author’s e-mail: tanabe@ikona.ed.shizuoka.ac.jp

５６７!2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.92, No.7 July 2016

ここでは並列化をコンパイラに指示するための規格…

OpenMP の概説をする．また，並列化を行う際に注意すべきこととして，スレッド間の負荷の均一化について触れる．

という処理をしている間に，影武者も全く同一の

２．
２．
２．
１ OpenMPc[0] = a[0] + b[0];C 言語や Fortran には，並列処理を行うための規格が存

c[1] = a[1] + b[1];在しない．そのため，並列処理機能を補う目的で１９９７年に…

業界標準規格として OpenMP が発表された．OpenMPの詳細を解説するには紙面が足りないので，ここではその代表という処理を行うことになり，並列化の意味がない．そこ的な使い方を示すに留める．詳しくはインテル㈱が提供しで，（３）の指示を書くことにより，本人は i＝０〜２
４の部ている「インテル C/C＋＋コンパイラー OpenMP 活用ガイ

分，一人目の影武者は i＝２５〜５０の部分，…，といった具合

ド」がコンパクトながら良くまとまっているので，参照しに負荷分散が行われる．この例からも，変数 i にプライベーてほしい［４］．

ト指定が必要な理由がよくわかるだろう．

OpenMP による C 言語でのプログラム例を示す．/* … */（５）はリダクションと呼ばれる演算のための書き方で，は，プログラム中での注釈を表す．

和を求める変数 sum にプライベート指定をしてしまうと，影武者たちが sum に求めた和の値を回収する機会を失/*_(1)_*/う．そこで，リダクション指定を行ったループでは，当該#pragma_omp_parallel_private(i)_/*_(2)_*/

ループ処理が終了したときに，（＋：sum）の後ろに指定｛されたプライベート変数＊２について（＋：sum）の前半に__#pragma_omp_for_/*_(3)_*/指定された演算（ここでは，"+"の加算）が行われる．他に__for_(i=0;_i<N;_i++)_{

も，リダクションでは乗算などを指定することもできる____c[i]_=_a[i]_+_b[i];が，可能な演算種には制限があるので，注意してほしい．

__}

最後に，
（６）に差し掛かると並列処理は終了し，影武者__#pragma_omp_for_/*_(4)_*/たちは消え去り，本人のみがその後の処理を継続する．以__for_(i=0;_i<N;_i++)_{

上が，OpenMP による並列化の基本的な考え方である．並____z[i]_=_x[i]_+_y[i];列化では，他にも，同期という重要な概念がある．これを__}きちんと管理しないと，走らせるたびに異なる結果が得ら__#pragma_omp_for_reduction_(+:sum)_/*_(5)_*/

れるといった問題が発生することになるので，文献などで__for_(i=0;_i<N;_i++)_{

チェックしてほしい［３］．

____sum_+=_c[i]_+_z[i];

２．
２．
２．
２ スレッド毎の負荷の均一化__}

例えば，以下のような４種類の処理があるとしよう．

}_/*_(6)_*/OpenMP では，並列処理をしない部分とする部分を行き来しながら処理を進める．まず初めに，プログラムは単一の処理の流れ（マスタースレッド）で処理が始まる（（１）の状態）．並列処理の開始を指示する宣言子に出会うと（（２）の（１）１〜１０の和を計算（２）１〜２０の和を計算（３）１〜３０の和を計算（４）１〜４０の和を計算

行），マスタースレッドである本人に加えて，それらの「影このような処理を，２分割して二つのプロセッサで並列に武者」が生成されて，本人と影武者による並列処理が始ま計算する場合，１番目のプロセッサに（１），（２），２番目のり，対応する閉じ括弧が現れる（６）まで続く．この影武者

プロセッサに（３），（４）と処理を割り振ると，２番目のプロを，スレーブスレッドと呼ぶ．
（２）の後半に指定されている

セッサの処理が終了する前に１番目のプロセッサの処理がprivate（i）は，並列処理中に含まれるループ処理でカウ終了し，１番目のプロセッサが遊休状態となる．すなわち，ンタとして使用される変数 i について，影武者用にも独立並列処理を行う場合には，プロセッサ毎の負荷を均一化しした変数 i を用意することを指示している．例えば，マスないといけないことがわかる．ベクトル化でもそうだが，

タースレッド一つとスレーブスレッド三つの計四つのス

計算のためのリソースをできるだけ多く活用するためにレッドで処理を行う場合，マスタースレッドの i に加えても，遊んでいるスレッドができるのは避けたい．これは，スレーブスレッド用に３個の i の格納場所が用意される．

コンパイラでは判断できないため，プログラマが明示的に

次に，（３）では i＝0，1，…，N−1 について行われるループ

OpenMP 宣言子 schedule を用いて指示する必要がある．通処理を，影武者と協力して分散処理することを指示する．

常 は， schedule（static） が 指 定 さ れ た も の と み な さ（３）の指示がないと，本人がれ，各スレッドにはあらかじめわかっているループの処理
回数を均等割した分の処理が割り当てられる．例えば，i＝c[0] = a[0] + b[0];０〜９９の処理を４スレッドで処理する場合には，１番目のc[1] = a[1] + b[1];＊２

スレッドには i＝０〜２
４，２番目のスレッドには i＝２
５〜４９

リダクション演算では，reduction 句の後ろの括弧内に指定された変数は，自動的にプライベートの扱いとなる．

５６８

Lecture Note

2. Let's Optimize a Program for Many-Core Accelerator

E. Tanabe et al.

という具合である．しかし，上記に示した例のように，各いるが，それよりも強い最適化が行われる．ベクトル化な

i 毎に処理の重さが異なっていると i での分割が平等でも実どに加えて，関数のインライン展開やループの展開（アン際の負荷は均等化されない．このような場合に scheduleロール）など，様々な最適化がコンパイラによって自動的（dynamic）を指定する．この場合には，１番目のスレッに行われる．-fp-model は，浮動小数点数の演算精度を制限

ドに i＝０，２番目のスレッドに i＝１，というように割りすることで，計算速度を向上させるオプションである．今振りながら，処理の終了したスレッドから順次，次の i が指回使用した-fp-model fast＝2 では，特に速度を優先した演

定され，処理を進める．staticの場合に比べて処理を配分す算が行われる．計算の精度は若干落ちるが，点渦系のようる手間はかかるが，遊休状態となるスレッドが減るため，な多体系では大きな問題が発生する可能性が低いことを経

並列化効率は結果として向上する．

験的に知っているため使用している．よって，計算モデル

２．
２．
３ FMA

によっては使用を控えた方が良い場合もある．-no-prec-div

FMA（Fused Multiply and Add）は，!"#!"という１は浮動小数点数の除算の精度を下げ，演算速度を上げるこ

次式を高速に計算する最適化技法である．上記のような積とを指示する．ただし，このオプションを指定すると，と和の計算が１回ずつ現れる式では，本来であれば加算命IEEE で決められた浮動小数点演算規則と異なる除算が行

令 と 乗 算 命 令 を１個 づ つ 使 用 し な け れ ば い け な い が，われるため，IEEE 準拠で除算を行った場合と計算結果がFMA では１命令で加算と乗算を行うことができる．FMA

異なることがある．よって，-fp-model と同様に計算結果にの適用の可否はコンパイラが自動的に判断してくれるが，影響を及ぼす可能性があるので，計算モデルによっては使複雑なコーディングでは見逃されることがあるので，やは用を控えた方が良い場合もある．なお，-fp-model と-no-

り１００％活用するには，わかりやすいプログラムとする必

prec-div は，-O3 が付いていなければ最適化が有効にならないため注意してほしい．

要がある．

２．
３．
２ 第２段階２．
３ 高速化の手順

ここでは，それぞれ

本節では，点渦系のシミュレーションコードを例にとって，具体的な高速化の手順（チューニング法）を紹介する．
チューニングは，３段階に分けて行う．第１段階は，研究者にとって負担の少ないコンパイルオプションのみによるチューニング，第２段階は，さらなる速度向上を狙ったプ
ログラムの書き換えによるチューニング，そして第３段階
は，第１，２段階を併用した最高速を狙うチューニングで（１）パディング
（２）配列のアライメント（３）pow 関数から２乗の掛け算への書き換え（４）２次元配列１本から１次元配列２本への書き換え（５）ベクトルループ内の配列要素の変数への書き換え（６）コア毎の計算負荷の均一化

ある．

に関する解説を行う．まず，キャッシュの利用率向上のた

ここで使用した機材／環境は，表１のとおりである．

めに，（１）と（２）を行う．パディングは配列要素へのアクセ

２．
３．
１ 第１段階スの際に発生するキャッシュの競合を緩和するために，本使用したコンパイラのオプションについて，概説する．

来必要とされる要素数に加えて，余分な要素を含めて配列-O0（オーゼロ）は，コンパイラによってデフォルトで行わ領域を確保する操作のことである．余分に確保する要素数

れる最適化を含めた全ての最適化を無効にするオプション

は，計算機や処理内容によって異なるので，いろいろと試である．普段使用することはないが，今回はコンパイラにす必要がある．今回，アクセラレータなしの場合は４要素，

よってどれだけ最適化がされているかを確認するために使アクセラレータ有りの場合は８要素を余分な要素として追用した．-no-vec は，ベクトル化のみを無効化するオプショ

加した．配列のアライメントは，配列を宣言する際に，ンである．-no-vec 以外にも個別に無効化できるオプション

All を配列要素数を指定する定数として，は多いが，アクセラレータの本質に関わりの高い最適化は_ _declspec(align(ALGN))

ベクトル化であるため，ベクトル化の有無の影響を調べ___________static_double_array[All];た．-O3 は，コンパイラによって最大限の最適化を行うことを指示する．コンパイラではデフォルトで‐O2 が付いてのように行う．この指定子を使うと，配列がプロセッサからアクセスしやすい境界に揃うように要素が並べられ，効

表１
プロセッサ
プロセッサ数プロセッサ当たりのコア数
主記憶容量
使用したアクセラレータ
OS
コンパイラ

率的に主記憶からキャッシュへデータが流れるようになる

環境．

ことが期待できる．アクセスしやすいメモリ境界は，プロIntel! XeonTM E5 2670
2
8
64GBIntel! Xeon PhiTM 7120A
CentOS release 6.7Intel! Parallel Studio XE 2015
Composer Edition for C++ Linux

セッサの設計特性に影響を受けるためマシンによって異なるが，アクセラレータの場合は６４バイト境界で始まるアド
レスにデータがあるとよい．
次に，（３）
について．関数 pow（a，b）は !" を計算する関数
である．かつては２乗の計算に pow 関数を使用すると高速になることが知られていた．そこで，Biot-Savart 積分に現れる２乗の計算に pow 関数を使用していたが，アクセラ
５６９

Journal of Plasma and Fusion Research Vol.92, No.7 July 2016

レータの場合には FMA との関連で，普通に a ＊ a と書い

しかし，この形式だとコア毎の負荷の分散が行いにくいよ

た方が速い場合が多いようである．よって，ここでは今まで，うので，次のように書き換えた．

pow（a，
２．
０）と書いていた部分を， a＊a と変更する．

for_(i=0;_i<N;_i++)_{（４）では，配列の書き方を変更する．例えば D を次元数，__for_(j=0;_j<i;_j++)

All を点渦数とすると，点渦の位置ベクトルを pos
［D］
［All］

____…

と宣言していたものを，posx
［All］
，posy
［All］
と書き換え

}

るよう に，２次 元 配 列 か ら 二 つ の１次 元 配 列 に 変 換 す

for_(i=0;_i<N;_i++)_{

る．２次元配列と１次元配列を比べると１次元配列の方が__for_(j=i+1;_j<N;_j++)

より少ない間接参照で当該要素にアクセスできる．またそ____…

れ以外にもコンパイラに「配列自身に依存性がない」と判}断をしてもらいやすいというメリットがある．
（５）は，ベクトル化が行われるループ内の配列要素を変これによって，コア毎の負荷分散がさらに促進され，結果＊数に書き換える操作である．例えば a［i］＝b［i］［
c i］と書としてスループットが向上する．

＊
いていたものを a＝b［i］［
c i］と書き換える．Intel コンパ

２．
４ ベンチマーク結果

イラでは，ベクトル化が行われるループ内にある変数は，ループ毎に（別々の変数として）プライベートに扱われる．

ここでは，前節で説明したチューニングを行った結果，

つまり，配列と同じように扱われるようになる．そのためどの程度高速化が達成できたのか定量的に示すため，ベン

本来なら配列を置かなければいけない部分を変数で置き換

チマーク結果を報告する．ベンチマーク共通の条件は，粒えることができる．変数と配列では，変数の方が間接参照

子 数 が１０３７０，ル ー プ の 回 数 が６０回 で あ る．最 終 的 に，の手間が一手間少ない＊３．すなわち，計算の途中結果など

Xeon E5 2670 ８コア×２での１６並列に対して約１０倍の高速

は配列ではなく，普通の変数に格納することによりメモリ化に成功した．これは，アクセラレータの理論ピーク性能

アクセスを高速化できる．

の５割程度に相当し，十分満足がいく結果といえる．

（６）で は，コ ア 毎 の 計 算 の 負 荷 を 均 一 化 す る た め に

２．
４．
１ アクセラレータを使用しない場合

OpenMPの宣言子で， schedule（dynamic） を用いる最適

ここでは，そもそも正しい挙動をするプログラムなのか化を行う．デフォルトの schedule（static）では，コア毎のを，アクセラレータを使用しないで確認する．初めに，プ

負荷などを考慮した並列化はできないが， schedule（dy-

ログラムにおいて，並列度を n 倍にしたときに，処理時間

namic）を明示的に指定することにより，実行時に動的に!
#になっているかを確認する．結果を，図１に示す．メが!

仕事の配分が行われるようになるため，並列化効率が向上

モリのアクセス競合等が発生するため，計算時間がきっちする． static と dynamic のどちらを使用するかは，動的り理論通りになることはない．それでも，スレッド数＝１に仕事の配分を行う手間と効率から判断すればよい．今回の場合と比べて，スレッ ド 数＝１６の 場 合 の 処 理 速 度 は

は，実行時に動的に並列化を行なうコスト以上に，スレッ14.5 倍となっており，今回作成したプログラムは，論理的ドへの負荷が均一化することによる速度向上の方が大きいに正しく並列化されていることがわかる．次に，スレッド

ため，dynamic を使用した．

数を１６に固定し，粒子数を変化させた場合の速度変化を確

２．
３．
３ 第３段階認する．結果を，図２に示す．縦軸が時間［秒］，横軸が点渦数［個］である．今回のホットスポット＊４である Biot-

第３段階では，第２段階で書き換えたコードに対して，第１段階で効果のあったコンパイラオプション-fp-modelと-no-prec-div を付けるチューニングを行い，どの程度高速化が可能かチェックを行った．さらに，Biot-Savart 積分の分母に表れる !!，!"という二つの位置ベクトルの差!!!!""がゼロとなる状態を取り除くため，チューニング
"前は，変数 i が固定された変数 j に関する内側ループで，j＝０から j＝i−1 まで と j＝i＋1 から j＝N−1 までの二つのループに分割していた．
for_(i=0;_i<N;_i++)_{__for_(j=0;_j<i;_j++)____…
__for_(j=i+1;_j<N;_j++)
____…

}

図１

＊３

配列は（レジスタではなく）主記憶に置かれるので，一般的に遅くなる．

＊４

プログラムで特に時間を要する重たい処理を，このように呼ぶ．

５７０

スレッド数をパラメータにした速度変化．

Lecture Note

2. Let's Optimize a Program for Many-Core Accelerator

E. Tanabe et al.

倍，-no-prec-div では 1.3 倍速くなる．
オプションはコンパイル時に指定するだけで良いので，研究者への負担が非常に少ないながら，オプション指定のみでアクセラレータなしでの１６並列と比べて 4.1 倍の速度を実現できた．
２．
４．
３ アクセラレータ＋プログラム書き換え
前節で説明したプログラムの書き換えを行った場合の結果を図４に示す．まず，（１）〜（５）の書き換えにより，書き換えを行う前で最も速い場合より 1.6 倍速くなっているのがわかる．さらに（６）の負荷の均一化を加えると，同様に
1.65 倍速くなる．速度向上効果は小さいように見えるが，粒子数が多くなるに従ってこの差は大きくなるので，書き換えを行う価値はあると考える．最終的に，（１）〜（６）の書
図２

き換えにコンパイラオプション‐fp-model fast＝２と‐no

粒子数をパラメータにした速度変化．点線は，N 2 に比例するシミュレーション時間の目安(傾き)を表す．

-prec-div を加えた場合には，（１）〜（５）の書き換えのみを行った場合と比較しても 1.5 倍速くなっている．

!

Savart 積分は，点渦数を ! とすると，計算時間は ! に比２．
４．
４ 最終結果

例する．点線が，!!に比例する時間計算量を表しており，

最後に，各段階でのベンチマーク結果の比較を行う．結

!

粒子数が多くなるにつれて ! に漸近していることから，

果を図５に示す．プログラムを書き換え，オプションを付

Biot-Savart 積分でのループの反復回数が間違っていない

けてチューニングを行った（４）は，アクセラレータなしのことを確認できる．なお，点渦数が少ない側で!!のスケー（１）に比べて 9.6 倍の速度に達する．プログラムの書き換え

ルからずれるのは，Biot-Savart 積分以外の部分の処理時間なしでオプションのみのチューニングを行った場合と比べに律速されていることを表している．

ても 2.3 倍となっている．手間を掛けず，オプションのみの２．
４．
２ アクセラレータ＋コンパイラオプション

チューニングでも，アクセラレータの有無で速度に 4.1 倍

通常のプロセッサ上で，期待通りの動作をしていることが確認できたので，ここから，アクセラレータを使用した
チューニングに移る．まず，研究者にとって負担の少ないオプションのみを付加するチューニングを行う．今回使用したオプションの組み合せは，（１）-O0（２）-O2-no-vec（３）-O3（４）-O3-fp-model fast＝2（５）-O3-no-prec-div（６）-O3-fp-model fast＝2-no-prec-divの六つである．結果を図３に示す．まず，コンパイラの恩恵を確認するために，（１）-O0 と（３）-O3 の比較を行う．両者 で は，約２６０倍 の 速 度 差 が あ る の が わ か る．す な わ

図３

コンパイルオプションをパラメータにした速度変化．

図４

プログラムの書き換えをパラメータにした速度変化．

ち，２６０倍に加速するだけの仕事を，研究者の代わりにコンパイラがしてくれたことになる．ただ，コンパイラも万能ではないため，研究者が最適化しやすいだろうと気を利かせてコードを書くと，コンパイラがこちらの意図を正確に判断できず，最適化をしてくれないことがあるため工夫をする時には注意が必要となる．
次に，（２）-no-vec と（３）-O3 を比較する．ベクトル化の有無で速度が６倍変化する．アクセラレータのベクトルレジ
スタは５１２ビットなので，完全なベクトル化を達成すると両者の差は８倍になるはずだが，残念ながら，そこまでの速度向上は望めなかった．
さらに，-O3 に加えて付加した（４）-fp-model fast=2 と（５）-no-prec-div の 結 果 を 見 て い く．-O3 と fp-model で は 1.4
５７１

Journal of Plasma and Fusion Research Vol.92, No.7 July 2016

たいのならば，コンパイラオプションだけを使うチューニ
ング方法もある．本章では，高速化に有効なコンパイラオ
プションとして-no-prec-div と-fp-model fast＝２の二つを紹介した．ただし，計算精度が落ちてしまうためモデルによっては導入を避けた方がよいかもしれない．コンパイラ
オプションだけでも速度が向上するが，さらなる高速化を期待するならば，プログラムを書き換える必要がある．まず配列の境界を６４バイトの倍数に揃えることで，メインメ
モリへのアクセス効率向上を狙う．次に，各種演算式を，できる限り，乗算と加算の組合せに書き換えて FMA が最大限利用できるようにする．また，１本の多次元配列を複図５数の１次元配列に書き換えると速くなることがある．ベク

最終結果．

トルループ内の変数は，ループ毎にプライベートに扱われるということを利用して配列を変数に書き換えてもメモリの開きが出る．そのため，オプションのみのチューニング

アクセスの効率を良くすることができる．最後に OpenMPか，プログラム自体の書き換えを行うかは，負担を考慮し宣言子で schedule（dynamic）を用いてコア毎の負担の均た上で選択すれば良いだろう．プログラムの書き換えは研一化を行う．コーディングによるチューニングとコンパイ究者にとって負担ではあるが，今回プログラムに施した

ラオプションによるチューニングは併用可能であるため，

チューニングは，ハードウェアやプログラミングへの深い

この両者を用いると更なる速度向上が期待できる．

知識を必要とせず，煩雑な手間が掛かるものではない．そのため，速度差が 2.3 倍になることを考え合わせると，プロ

参考文献

グラム自体を書き換えるチューニングは施す価値があると［１］J. Jeffers, J. Reinders，（訳）
すがわらきよふみ，エクセル
ソフト株式会社：インテル Xeon Phi コプロセッサーハ
イパフォーマンス・プログラミング（カットシス テ
ム，東京，２０
１
４）．
［２］D.A. Patterson，（訳）
成田光彰：コンピュータの構成と設計第５版（上下）（日経 BP 社，東京，２０１４）．
［３］北山洋幸：OpenMP 入門（秀和システム，東京，２００９）．
［４］インテル ! ： http : / / jp. xlsoft. com / documents / intel /compiler/526J-001.pdf.

考えている．

２．
５ 結論
アクセラレータの高速化技法には，ベクトル化，並列化，
FMA 等がある．アクセラレータをより高速化させるためには，これらを最大限運用できるようにコーディングを工夫する必要がある．もし，コーディングの負担を小さくし

５７２

J. Plasma Fusion Res. Vol.92, No.8 (2016)６
１２‐６２２

講座 今日からはじめるメニーコアアクセラレータ

３．IFERC-CSC におけるメニーコアアクセラレータとその利用例
3. The Many-Core Accelerator and Examples of Its Use at IFERC−CSC
中 島 徳 嘉，宮 戸 直 亮１），内 藤 裕 志２），佐 竹 信 一３），保 坂 和 樹３），功 刀 資 彰４）
NAKAJIMA Noriyoshi, MIYATO Naoaki１）, NAITOU Hiroshi２）, SATAKE Shin-ichi３）,
HOSAKA Kazuki３）and KUNUGI Tomoaki４）
核融合科学研究所，１）量子科学技術研究開発機構，２）山口大学大学院創成科学研究科，
３）

東京理科大学基礎工学部，４）京都大学大学院工学研究科（原稿受付：２
０
１
６年５月２
０日）

最近の大型並列計算機では，性能向上および省電力化のため，一つのチップに多数の演算コアを搭載する傾向にある．青森県六ヶ所村にある国際核融合エネルギー研究センター（International Fusion Energy ResearchCentre; IFERC）の３副事業の一つである計算機シミュレーションセンター（Computational Simulation Centre;
CSC）の大型並列計算機 Helios では，CPU のみの本体システムに加えて，多数の演算コアを持つコプロセッサを搭載する増強システムを運用している．本章では，この増強システムの利用環境，増強システムを実際に利用した２次元 PIC シミュレーションおよび MHD 乱流の直接数値計算について紹介する．
Keywords:
BA agreement, IFERC, CSC, simulation project, coarse-grained parallelism, fine-grained parallelism,
streaming algorithm, direct numerical simulation, MPI+OpenMP, MIC

３．
１ IFERC-CSC の MIC システムの環境紹介

速・並列なデータの入出力が可能となっている．

国際核融合エネルギー研究センター（International Fu-

MIC シ ス テ ム で は プ ロ グ ラ ミ ン グ ツ ー ル と し て，sion Energy Research Centre; IFERC）の計算機シミュレー

OpenMP（Open Multi-Processing）4.0 を含む Intel コンパ

ションセンター（Computational Simulation Centre; CSC）
は

イラ（C/C＋＋，FORTRAN），Allinea DDT，Intel VTune

日 欧 間 の よ り 広 範 な 取 組 を 通 じ た 活 動（Broader Ap-

といったデバッガが利用可能である．また，Intel MKL，

proach; BA）協定に基づいて実施されている事業で，青森

FFTW，PETSc，HDF5，MAGMA MIC，NAG，NetCDF

県六ヶ所村に大型並列計算機 Helios（愛称：六ちゃん）がという数値計算でよく使われるライブラリ類も MIC シス

設置され運用されている．この Helios 計算機は本体システ

テム用のものが用意されている．Message Passing Inter-

ム（４５００ノード）と増強システム（１８０ノード）から成り，

face（MPI）ライブラリとしては本体システムと同様に，増強システムに Intel 社のメニーコア（Many Integrated

IntelMPI と BullxMPI の２つが利用可能となっている．

Core;MIC）アーキテクチャをベースにした Xeon PhiTM コ

ユーザーは，図２に模式的に示されているように，（１）プロセッサが搭載されている．このため増強システムは
MIC システムとか MIC ノードと呼ばれている．
Helios 計算機の導入の経緯などは２０１５年１１月号掲載の解説記事［１］に詳しいのでそちらを参照されたい．ここでは
MIC システムの環境について簡単に説明する．MIC システ
ム の１ノ ー ド は，２つ の Intel 社 製 CPU（Xeon，SandyBridge EP，８コア，2.1 GHz）および２つの Xeon Phi コプ
ロセッサ（Xeon Phi 5110P，６０コア，1.05 GHz）と，合計で６４GB のメモリを搭載している（図１）．MIC システム全体
としての理論ピーク性能は４１２TFLOPS で，linpack 性能として 225.1 TFLOPS が得られているが，これは２０１５年１１月の TOP500 リスト［２］で４５６位に位置する．この MIC シス
図１

テムは，本体システムとファイルシステムを共有し，高１つの MIC ノードの構成．

corresponding author’s e-mail: nakajima @nifs.ac.jp

６１２!2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

3. The Many-Core Accelerator and Examples of Its Use at IFERC-CSC

N. Nakajima et al.

Phi のプログラミング方法についての一般的な解説書としては文献［３］があり，本章でも参考にした．
３．
２．
１ Xeon Phi でのスレッド並列化の指針
科学技術計算の高速化には，並列化が必須である．並列化には，プロセス並列とスレッド並列がある．プロセス並
列は分散メモリモデルに対応しており，各プロセスは独立なメモリを持ち，独立に動くが，必要に応じて，異なるプ
ロセスの持つメモリ情報を交換（通信）したり，同期を取ったりする．通常は MPI を用いてプログラムする．ス
レッド並列は共有メモリモデルに対応しており，各スレッドから共通のメモリ空間への読み書きができる．通常マス
タースレッドが全体を統括し，必要に応じて，自分を含む複数の管理下のスレッドに実行命令を送る．並列化には，
自動並列化か OpenMP を用いる．複数のコンピュータ環境を利用する場合は，可搬性に優れた OpenMP の利用が推奨される．
単一のパーソナルコンピュータ（PC）のユーザーは，マルチコアの性能を引き出すために，プロセス並列かスレッ
図２

ド並列かのどちらか一つを採用する．プロセス並列を利用

コプロセッサを利用する３つのモード．

する場合，共有メモリは仮想的に分散メモリとして利用される．PC クラスタやスーパーコンピュータのユーザーは，

オフロードモード，（２）ネイティブモード，（３）シンメト

純粋なプロセス並列モデルを用いるか，プロセス並列とスリックモードの３つの形態で MIC システムのコプロセッ

レッド並列のハイブリッドモデルを用いる．スレッド並列

サを計算に用いることができる．

部分で計算速度の十分な高速化を得るためには，スレッド（１）オフロードモード

間の独立性の極めて高い（場合によっては新しい）アルゴ

MIC ノードの中のホスト CPU でプロセスを実行し，リズムを採用することが必須である．

コードの中で指示行により明示した部分だけ，コプロ

アクセラレータとして Graphics Processing Unit（GPU）

セッサを利用するモード．

や Xeon Phi を用いる場合は，スレッド並列が唯一の選択肢（２）ネイティブモード

である．プログラム開発をできるだけ効率良くするためにコンパイラでコプロセッサ用のバイナリを作成し，コは，マルチコアシステムで並列性の高い（コア数に比例しプロセッサのみで計算を行うモード．

て高速化される）アルゴリズムを採用したプログラムを（３）シンメトリックモード

OpenMP で書き，デバッグを済ませておくことが重要であMPI によってホスト CPU，コプロセッサで複数のプロ

る．こうしたプログラムを標準プログラムにして，GPU

セスを実行し，計算を行うモード．このモードでは異の場合は，例えば CUDA FORTRAN を用いて書き換え，なるノード上の MIC（Xeon Phi）同士でもホストを介

Xeon Phi の場合はそのまま利用（ネイティブモードの場

さずに MPI を使った直接通信が可能．

合）するか，指示行を追加することにより書き換える（オこのように MIC システムの利用法には多様性があり，そフロードモードの場合）．オフロードモードの場合は，６０

れに応じて，プログラム作成の難易度やコードの実行速度

コアの内１コアはホストとの通信に使用されるため，使用も変化することになる．ユーザーの方々が今後 MIC システ

できるコアは５９コアになる．１コアは４ハードウェア・ス

ムを活用する際の具体的参考例として，３．
２節では主としレッドに対応するため，１個の Xeon Phi で，ネイティブてネイティブモードを用いた２次元 PIC シミュレーション

モードの場合は２４０スレッドまで，オフロードモードの場を内藤が，３．
３節ではシンメトリックモードを用いた MHD

合は２３６スレッドまでの計算が可能になる．

乱流の直接数値計算を佐竹，保坂，功刀が紹介する．これ

ハイブリッド並列の場合は，MPI を併用する．ネイティ

らの具体例を参考として，MIC システムの効率的な利用がブモードの場合は，シンメトリックモードを用いる．この広まることを期待している．

場合６０コアの内，１コアはホストとの通信に使用されるため，スレッド並列で使用できるのは５９コアである．オフ

３．
２ ２次元 PIC シミュレーション

ロードモードの場合は，そのまま MPI を併用することがで

本節では，２次元 PIC（Particle-In-Cell）コードの並列化

きる．

を例に Xeon Phi 対応プログラム作成の指針と応用例につGPU をアクセラレータとして使用するのと比較して，いて解説する．説明は，IFERC-CSC の MIC システムで使

Xeon Phiを利用する最大の利点はマルチプロセッサ対応の用されている Xeon Phi 5110P を具体的な対象としている

プログラムがあり，それが十分な並列化性能を持っていれが，その他の Xeon Phi についても本質は同じである．Xeon

ば，ほとんど無修正で Xeon Phi で実行できることである．
６１３

Journal of Plasma and Fusion Research Vol.92, No.8 August 2016

ただし，プログラムの一部に並列化が不可能もしくは並列

従 来，PIC コ ー ド の 並 列 化 に は，粗 粒 度 の 並 列 性

性に著しく劣る部分があって，全体を並列実行した場合に（coarse-grained parallelism）が用いられてきた．これはプ

その部分の実行時間が無視できない場合は，オフロード

ロセス並列に対応している．場の量に複数のコピーを使用

モードを用いる．この場合，ホストプロセッサと Xeon Phi

する粒子分割や，場の量を複数の領域に分割する領域分割

コプロセッサでそれぞれ得意な分野の計算を分担することを用い，主に MPI を用いてプログラムされる．詳しくは本になる．

誌の講座［６］を参照されたい．先に示したように，従来の３．
２．
２ Xeon Phi での PIC コードの並列化

PIC コードの手法ではスレッド並列の部分の並列化性能を

ここでは，２次元 PIC コードを例にネイティブモードで引き出すのは困難である．筆者の経験では，すべてプロセの Xeon Phi の利用について述べる．PIC コードは，荷電粒ス並列でするよりは，２スレッドか４スレッドのスレッド

子の運動を，荷電粒子自身が作る電磁場（もちろん外部磁並列を併用したハイブリッド並列の方がわずかに高速で場や外部電場を含めて良い）による電磁力のもとで追跡すあった．

ることにより，プラズマのシミュレーションを行う．第一マルチスレッドやアクセラレータに対応して，PIC コー

原理シミュレーションであるため，コンピュータに対するドの性能を引き出す手法の一つに，細粒度の並列性（fine-

負荷は巨大であるが，最新のコンピュータ環境を利用する

grained parallelism）を利用する方法がある．ここで解説すことにより，核融合プラズマや実験室プラズマの物理現象

る手法は，Decyk と Singh による論文［７，８］に詳しく記述さを解明したり予測したりする重要な武器になっている．

れている．これらの論文では，GPU での高速化がメイン

PIC コードによるシミュレーションの教科書としては文献

テーマであるが，マルチコアやメニーコアへの拡張性につ［４］がある．また，本誌にも講座［５］がある．

いても議論されている．

PIC コードでは，粒子と場の量を取り扱う．粒子とはプ

PIC コードでは，任意の粒子は，その粒子の存在するセラズマを構成する電子やイオンの荷電粒子を示し，シミュ

ルのみと相互作用する．このため，粒子全体を一つの配列

レーション空間を自由に動きまわる．場の量は電場と磁場で表すのではなく，セルごとの配列を用意する．これは粒を示し，空間メッシュの格子点でのみ計算される．ここで子がセルごとにソートされていることに対応している．まは，直角座標で，各方向で等間隔にメッシュが構成されてた，全体の場の量を表すメッシュ配列に加え，セルごとのいるとする．空間メッシュを構成する最小単位をセルと呼ぶ．

配列を用意する．SOURCE と PUSH の計算では，従来の任意の時間で任意の荷電粒子を取り出すと，空間メッ

PIC コードが粒子番号に対応する DO ループで表されたの

シュ上のどこかのセル内に存在する．この荷電粒子が感じと比較して，新しいアルゴリズムではセル番号に DO ルー

る場の量は，近傍の格子点の場の量からの内挿で計算されプが対応する．また内側の DO ループとしてセル内の粒子

る．線形補間を用いる場合はセルの頂点の場の量のみが使に対するループを追加する．各セル内の粒子の総数を与え用される．以下では線形補間の場合に限定して解説する．

る配列も用意する必要がある．計算はセルごとに独立であ粒子の位置と速度は，ニュートン・ローレンツの運動方程るため並列性は極めて高い．GPU の場合は，１スレッドが式を差分化したものを用いて１時間ステップ進める．この１セルに対応している．Xeon Phi の場合は，複数のセルを過程を PUSH と呼ぶことにする．PUSH では，粒子番号の１スレッドで対応することになるが，セル単位で独立性が順に処理を進めるが，各粒子が引用する場の量へのアクセあるため，スレッドごとの独立性も保証されている．また，

スはメモリ空間でランダムになる．このことは，キャッメモリーアクセスの場合の局所性も保証されるのでキャッシュの有効利用という点で障害になり，コードの高速化性

シュの性能も引き出しやすい．

能の劣化の原因になる．また，この部分を並列化した場合，新しいアルゴリズムを文献［７，８］ではストリーミング・

メモリへの同時アクセスが生じやすく，速度低下の原因にアルゴリズムと呼んでいる．PUSH の処理を終了すると，なる．

セルの境界を超える粒子を，新しく粒子の所属するセルへ場の量を計算するためには，格子点上の電荷密度と電流移動させる過程（SORT）が追加される．１時間ステップ

密度を計算する必要がある．各粒子の電荷と電流を，セルに，複数のセルを横切る粒子はないと仮定して問題ない．

の頂点に対応する格子点に分配する．分配公式としては，

以 上 の 並 列 化 を セ ル 並 列 と 呼 ぶ．セ ル 並 列 の 場 合，

逆線形補間を用いる．この過程を SOURCE と呼ぶことに

SOURCE と PUSH の並列化性能は優れているが，追加されする．SOURCE では，場の量にランダムにアクセスし，た SORT の計算時間が無視できなくなる．この問題を解決

データを書き換える（加算する）必要があり，高速化の妨するために複数のセルをまとめてタイルとしたタイル並列げになっている．また，並列化した場合は，同時に同一のを利用する．

データを書き換える場合の結果は保証されていない．ハー２次元でのセルとタイルの例を図３に示す．タイルの例

ドウェア的またはソフトウェア的に並列化することは可能

では "方向４（!"!4），#方向４（!#!4）の場合を示しであるが，速度低下の原因になっている．

ている．タイルを利用するとセルの場合と比較して，全体電荷密度と電流密度から場の量を計算する過程を FIELDとして，領域境界を横切る粒子の数が減るため計算時間がと呼ぶことにする．PIC コードは，SOURCE-FIELD-PUSH

短縮される（タイルを使用すると，タイル境界を除くタイを順に計算することにより１時間ステップ進めることがでル内のメッシュを横切る粒子に対する SORT の処理が必要

きる．

なくなると考えれば理解しやすい）．OpenMP で書かれた
６１４

Lecture Note

3. The Many-Core Accelerator and Examples of Its Use at IFERC-CSC

図３

N. Nakajima et al.

２次元でのセルとタイルの例．

SORT のプログラム例が文献［８］に記載されている．タイ
ルを利用した場合，タイルサイズが大きくなるにつれて，
SOURCE と PUSH の計算時間は増大する傾向にある．この

図４

マルチコアに対する２次元 PIC コードの並列化性能．

図５

マルチコアに対する２次元 PIC コードの高速化因子．

ため，タイルサイズには最適値がある．
３．
２．
３ Xeon Phi での２次元 PIC コードの実行例
ここでは，２次元静電近似コードの Xeon Phi での実行例に つ い て 記 す．プ ロ グ ラ ム は FORTRAN で 書 か れ，OpenMP 4.0 により並列化されている．場の量は，
２５６×２５６のメッシュで表される．各セル当たりの粒子数は平均して１００個，全粒子数は６，５５３，６００である．#と $の両方向には周期境界条件を用いる．外部磁場はない．電子のみを考え，イオンは電荷を中和するための一様なバックグラウンドとする．粒子に関する量は，位置 #，$と速度 "#，"$がある．
場の量としては，電場の #，$成分，電荷密度と静電ポテン
シャルがある．粒子の位置と速度を１ステップ進めるためにリープ・フロッグ法が用いられている．電子の空間分布は一様で，速度分布はマクスウェルである．空間はメッシュ幅で，時間はプラズマ振動数の逆数で規格化されている．規格化されたデバイ長（規格化された電子の熱速度と同じ）は１であり，規格化された時間ステップ幅は 0.1 である．１０００ステップの計算時間を比較する．計算はすべて倍

短縮はわずか（約８パーセント）である．使用した各スレッ精度で実行した．結果の一部は，本学会誌に掲載されたプドで，PUSH の計算時間の割合が最大で，SOURCE，SORT

ロジェクトレビュー［９］の１８９‐１９０ページに示されたものとが続く．FIELD の計算時間の割合は無視できる程度であ

ほぼ同様であるが，本小特集の執筆に際して高速フーリエる．図５は，SOURCE，PUSH，SORT，FIELD に対して，変換（Fast Fourier Transform; FFT）の並列化を加え，再それぞれスレッド数１の場合と比較した高速化因子を示し計算した結果を示す．またホストプロセッサでハイパースている．８コア（１６コア）で比較するとSOURCEで7.4（7.1）

レッディングを利用した計算を追加した．

倍，PUSHで7.4（8.1）
倍になっているので並列化性能が高い

図４と図５は，IFERC の Helios 増強システムの単一のホことが実証されている．SORT と FIELD についても 5.8

ストプロセッサ（Xeon，８コア）での２次元 PIC コードの

（7.8）
倍，6.4
（8.0）倍になっている．全体としては 7.1（7.7）倍並列化性能を示している．タイルサイズは，Xeon Phi の場の高速化が達成されている．SOURCE の場合のみハイパー

合に最適になる !#!!，!$!!を用いた．図４は，使用しスレッディングを適用して少し遅くなっているのは気になた ス レ ッ ド 数 に 対 し て の SOURCE，PUSH，SORT，

るところである．

FIELD の部分の計算時間を積み上げ棒グラフで示していマルチコアで高い並列性能が証明されたので安心してる．８スレッドまでは物理コアに対する計算でありスレッ

Xeon Phi に適用できる．図６と図７は，Xeon Phi での並列ド数の増加に対応して高速化されている．１６スレッドの場化性能を示している．図６は，使用したスレッド数に対す合は１物理コアを２論理コアとみなすハイパースレッディ

る計算時間の積み上げ棒グラフを示している．スレッド数

ングを用いているが，８コアの場合と比較して計算時間のの増加に従って，順調に高速化されていることがわかる．
６１５

Journal of Plasma and Fusion Research Vol.92, No.8 August 2016

図８

メニーコア（２４０スレッド）とマルチコア（１６スレッド）の２次元 PIC コードの実行時間の比較．

Xeon Phi の計算時間の高速化因子を示している．SOURCE
図６

メニーコアに対する２次元 PIC コードの並列化性能．

で 2.0 倍，PUSH で 2.2 倍，SORT で 2.2 倍に高速化されている．FIELD は 0.28 倍と遅くなっているが全計算に対する比率が小さいため問題ではない．全体で 1.9 倍の高速化が達成されている．これは，１スレッドで比較すると Xeon
Phi の方が 5.8 倍遅いが，スレッド数の差とハードウェア・
スレッドによって逆転して高い高速化が達成されているこ
とになる．
３．
２．
４ まとめ
２次元静電近似 PIC コードを例に，Xeon Phi での高速化手法について解説した．実行時間を比較すると，Xeon Phi
を用いると，ホストプロセッサに対して約２倍の高速化が達成されることを実証した．筆者は同じ計算を IFERC-CSCの 新 し い GPU 並 列 シ ス テ ム の 構 成 要 素 で あ る GK210（TESLA K80 に２個搭載されている）と核融合科学研究所の新しいプラズマシミュレータの構成要素である
４９６コアを持ち最SPARC６
４TMXIfx で実行した．GK210 は２大理論性能は0.935 TFLOPSである．SPARC64 XIfxはマルチコアであり，３
２コアに加えて２個の補助コアを持ち，

図７0.908 TFLOPS で あ る．Xeon Phi 5110P の 理 論 性 能 はメニーコアに対する２次元 PIC コードの高速化因子．

1.011 TFLOPS である．Xeon Phi 5110P と GK210 はアクセ
図７はスレッド数の増加に対する高速化因子を示していラレータであり，SPARC64 XIfx はホストプロセッサである．６０スレッドで比較すると，SOURCE で 53.3 倍，PUSH

る違いがあるが，理論性能はほとんど同じである．またで 55.8 倍と高い並列化性能を示している．SORT では 45.3

SPARC64 XIfx は OpenMP で書かれたプログラムがそのま倍である．FIELD は 22.4 倍と並列化性能は少し劣るが，全ま利用できる．結果として，SOURCE，PUSH，SORT，体の計算時間に対する割合が小さいため問題ではない．ま

FIELD の計算時間の割合はそれぞれ異なった傾向を示し

た，FIELD では FFT の計算が支配的であり，３
２スレッド

たが，全体の実行時間はほとんど同一であった．メニーコ

以上で高速化は飽和している．プログラム全体では 51.8

アに対してマルチコアのコア数も増加傾向にあるため，次倍に高速化されている．

世代で差がどのように変化するか興味があるところであ１２０スレッドと２４０スレッドの計算は，ハードウェア・ス

る．また GPU についてはプログラムの難しさはあるが，すレッドが用いられているためか，マルチスレッドの場合とでにプログラムを開発済みのプログラマにとっては問題な比較して顕著な高速化が測定されている．２
４０スレッドでい．それぞれどのような進化をするのか楽しみである．

比較すると SOURCE が 76.8 倍，PUSH が 110.1 倍，SORTが 84.0 倍と高速化されている．FIELD は 19.6 倍となっているが，全体の計算時間に対する割合が小さいため問題ではない．全体では 86.2 倍の高速化を示している．
図８は，１６スレッ ド の ホ ス ト プ ロ セ ッ サ と 比 較 し た
６１６

Lecture Note

3. The Many-Core Accelerator and Examples of Its Use at IFERC-CSC

N. Nakajima et al.

３．
３ マルチ MIC＋マルチ CPU で構成されたヘテ
ロジニアス計算機を使った MHD 乱流の直接
数値計算

non-slip 境界条件を，流れ方向とスパン方向には周期境界

３．
３．
１ はじめにド法である．計算格子は，壁垂直方向速度のみにスタガー条件を課した．空間の離散化には，２方向にスペクトル法
を用いて壁垂直方向に差分法を用いるいわゆるハイブリッ

近年のスーパーコンピュータは GPU や MIC のようなアド格子を用いる．粘性項の時間積分には，陰解法である

クセラレータが使われているものが多数を占めている．こCrank-Nicolson 法を用い，他の項には陽解法である三次精

れらは MPI を用いた従来の並列技術に加え，各コアでス

度 modified Runge-Kutta 法を用いる．

レッド並列を行っている．オフロードモードではGPUを用

３．
３．
２．
１ 並列化の実装いた並列計算をアクセラレータ（コプロセッサ）に処理を

粘性項のHelmholtzの式と圧力Poisson方程式は，フーリさせることが必要となる．しかし，MIC と呼ばれる Intel

エスペースの三重対角行列解法を用いて解いている．図９

Xeon PhiコプロセッサはGPUとは異なり，コードを書き換で示すように並列アルゴリズムとして !
-"面内で FFT をかえることなく実行できる．このモードでは OS が CPU 側のけ，壁垂直方向に三重対角行列解法による解法を必要とす

ホストプロセスで実行されるため，コプロセッサ側で処理

る．Poisson 方程式は図９のような転置を必要とし，メモリする必要がない．したがって，MPI プロセスを直接コプロ

パーティションが異なる２つのメモリ空間を用いる．毎ス

セッサで実行することができるようになる．その場合，コテップで対流項計算時（パディング法）に転置が必要となプロセッサに MPI のランクi が設定でき，ピア・ツー・ピ

る．等間隔計算格子と周期的境界条件は主流方向で使われア通信を行うため，独立した計算ノードのように働く．しる の で，フ ー リ エ 変 換 を 用 い る こ と に よ っ て３次 元

たがって，プロセッサとコプロセッサで構成されるヘテロ

Poisson 方程式を１次元の式として扱うことができる．行

ジニアスシステムが可能となる．つまり，MPI+OpenMP列の計算において，三重対角行列解法は，波数ごとに使用で記述されたプログラムコードを変更することなく，コプしている．

ロセッサの有無にかかわらず MPI アプリケーションを使

３．
３．
２．
２ プログラミングの方針用が可能である．本節で紹介する乱流直接数値シミュレー

MPI+OpenMPを適用したFORTRANで書かれたプログ

ションのプログラミングコードはMPI + OpenMPによる並ラムコードをネイティブモード（MIC）とヘテロジニアス

列化をしている．その際の効率的なヘテロジニアスシステ

計算（CPU+MIC）で使う．本プログラムコードは FOR-ムのためのコード設計，例えばサブミッションシェルにおTRAN９０で書かれており，MPI ライブラリを参照していけるノード数，CPU 数，および MIC 数の組み合わせと計算

る．３次元配列の計算に対し，外側の DO ループを MPI

速度について報告する．
３．
３．
２ DNS コード
すでに我々は，CPU+MIC を用いて円管内乱流 DNS コー
ドの 並 列 化 に 成 功 し て い る［１０‐１２］．本 節 で は，Satake
et al.［１３］により開発されたチャネル内乱流 DNS コードを改良したものを用いて，Helios 計算機を使った速度向上についてベンチマークを行った．まず用いたチャネル内乱流
DNS コードの計算手法について説明する．３次元乱流の運動 を 規 定 す る 支 配 方 程 式 は，直 交 座 標 系 で の Navier-

図９

Stokes 方程式，熱伝導方程式，連続の式である．壁面には計算領域に対するメモリーの分割方法と転送．

図１
０ シンメトリックモードにおけるランクの構成．

i ランクとは，MPI でのプログラムの実行単位であるプロセスの識別番号で，０から順に連続した整数が割り当てられる．

６１７

Journal of Plasma and Fusion Research Vol.92, No.8 August 2016

による並列化を行い，その内側の２つのループに関してで使用しているため，通常は使用することができない．

OpenMP による並列化を行った．このルールは FFT や軸よって MIC の最大スレッド数は２３６スレッドとなる．

交換ルーチン等，全てに適用している．

３．
３．
２．
４ ヘテロジニアス計算を行うための実行シェル

３．
３．
２．
３ Helios 計算機のヘテロジニアスクラスタリング

ホスト CPU と MIC を同時に使うヘテロジニアス計算をシステム行うには使用するノードと CPU 及び MIC のプロセッサを

シンメトリックモードは，図１０のようなヘテロジニアス

シェルで指定する必要がある．１ノード内で 2CPU，2MIC計算である．

を稼働させ，１つのプロセッサ内で２つランクを立てた時の例を図１１に示す．また，ここでは主要な記述以外は省略

このモードではアプリケ ー シ ョ ン プ ロ セ ス が ホスト
CPUとXeon Phiコプロセッサの両方で実行される．それらしている．

は通常，MPI を介して通信する．この実行環境はヘテロジ

まず NODELIST で使用するノードを決定する．この記

ニアスクラスタ環境において同一クラスタ内でも別ノー述では１６個のノードを使用するため，-I0-15 と記述していド と し て 扱 う こ と が で き る．ス レ ッ ド 定 義 はる．次に for ループで使用するノードと CPU 及び MIC を

KMP̲AFFINITY によって指定する．スレッドは０からMachinefile に書き込む．CPU を使う場合はノードの番号

２３９の範囲で KMP̲AFFINITY によって書くことができのみ書き込み，MIC を使用する場合はノード番号のあとにる．し か し，０，２３７，２３８，お よ び２
３９の ス レ ッ ド は OS

-mic0 もしくは-mic1と書く．今回のMachinefileへの書き込

図１
１ 実行シェル（hybrid2.sh）．

６１８

Lecture Note

3. The Many-Core Accelerator and Examples of Its Use at IFERC-CSC

N. Nakajima et al.

みは for ループで１６回書くことを４パターン，２回繰り返のみをif文で判定している．ランクが１
６より小さい場合は，し行っている．つまり使用するノードが１６個で，１つのKMP̲AFFINITY でスレッド０から３を指定し，mic̲tmp

ノード内で CPU を４個，MIC を４個使うことになる．よっ内の pipe.out をホスト CPU で実行する．ランクが３
２より小て，ランクが１２８個設定される．なお，ランクは Machinefile

さい場合は，スレッドを KMP̲AFFINITY で指定し，mic̲

に書かれた順番通りに振られていく．次に CPU 用の実行

tmp 内の pipe.mic を MIC で実行する．ランクが４８より小さファイル pipe.out と MIC 用の実行ファイル pipe.mic，実行い場合は，KMP̲AFFINITY でスレッド４から７を指定

用のシェルの pipe.hyb2.sh を MIC 用のディレクトリである

し，mic̲tmp 内の pipe.out をホスト CPU で実行する．ラン

mic̲tmp にコピーする必要がある．そのため，clush というクが６４より小さい場合は，スレッドを KMP̲AFFINITY

コマンドを使いコピー先のノードを指定し，各ノードにで指定し，mic̲tmp 内の pipe.mic を MIC で実行する．同様

mic̲tmp を 作 り，そ の 中 に 上 記 の フ ァ イ ル を コ ピ ー すにして，１２８ランクまで指定する．なお，１つのプロセッサ

る．ジョブ実行コマンドは mpiexec.hydra である．実行で２つのランクを立てるときには使用するスレッドが重複

ファイルは直接実行するのではなく，別のシェルで実行すしないように注意しなければいけない．さらに MIC の場合る必要がある．この時のシェルの記述内容を図１２に示す．

は OS で使用するスレッドも使ってはいけない．今回の場

また，ここでは主要な記述以外は省略している．

合は mic0，mic1 共に０，１１９，１２０，２３９が OS で使用されこのシェルは Machinefile の記述によってどのプロセッ

ている．
３．
３．
３ 結果と考察

サでどのファイルを実行するか判定している．Machinefile

表１にある２種類のレイノルズ数を使い，各レイノルズに mic と書いているかの判定はしておらず，ランクの番号

図１
２ 実行シェル２（pipe.hyb2.sh）．

６１９

Journal of Plasma and Fusion Research Vol.92, No.8 August 2016

数に応じたメッシュ解像度でそれぞれのベンチマークをた．４ノードでは CPU+MIC がホスト CPU の速度を超え行った．

た．図１６に１０２４×１０２４×７６８において合計ランク数を１２８に

これらのデータは Satake et al.［１３］に基づいている．プ

固定した時，ホスト CPU，MIC，CPU+MIC を様々な組み

ログラムのベンチマークは，経過時間［sec］を倍精度実数

合わせで稼働させた時の速度比較を示す．表３は図１６におで返す関数である mpi̲wtime（ ）を用いて計測した．

ける１つのCPU，１つのMICあたりのランク数の組み合わ

３．
３．
３．
１ ヘテロジニアス計算（ホスト CPU とマルチのせを示している．

MIC）の評価

ランクの数の合計はどのパターンにおいても１２８となる．

比較を行うためにノード数を同じにした CPU と MIC

全てのケースにおいて MIC 内のスレッド数は１２０に固定しの結果も示した．図１３は，１２８×１２８×１２８での CPU+MIC

た．３２ノードでは CPU+MIC がホスト CPU の速度を超えの ノ ー ド 数 に 対 す る 速 度 向 上 率 を 示 し て い る．図１４は，１０２４×１０２４×７６８での CPU+MIC のノード数に対する速度向上率を示している．CPU のみで行った計算が最速
で，次は CPU+MIC，最も遅い結果が MIC のみとなった．
図１５にランクの数を固定した時，１２８×１２８×１２８においてホスト CPU，MIC，CPU+MIC を様々な組み合わせで稼働させた時の速度比較を示す．表２は図１５における１つのCPU，１つのMICあたりのランク数の組み合わせを示している．
ランクの数の合計はどのパターンにおいても３２となる．
全てのケースにおいて MIC 内のスレッド数は１２０に固定し２８
図１
５ ３
２ランクの場合のノード数変化による計算時間（１
２
８×１
２
８）
．
×１
表２

128×128×128 における CPU と MIC の組み合わせ．

Nodesdata/CPUs
data/MIC0
data/MIC1

３
２
１
０
０

CPUs
１
６ ８
２ ４
０ ０
０ ０

４
８
０
０

３
２
０
１
０

MICs
１
６ ８
０ ０
１ ２
１ ２

CPUs+MICs
１
６ ８ ４
１ ２ ４
１ １ ２
０ １ ２

図１
３ ノード数変化による計算時間（１
２
８×１
２
８×１
２
８）
．
表１

レイノルズ数とメッシュ数．

Reynolds number
150
1100

number of mesh points (x y z)
128x128x128
1024x1024x768

図１
６ １
２
８ランクの場合のノード数変化による計算時間（１
０
２
４×
１
０
２
４×７
６
８）
．
表３

1024×1024×768（128 ranks）における CPU と MIC の組み
合わせ．
CPUs

Nodes
data/CPUs
data/MIC0
data/MIC1

図１
４ ノード数変化による計算時間（１
０
２
４×１
０
２
４×７
６
８）
．

６２０

１
２
８
１
０
０

６
４
２
０
０

３
２
４
０
０

MICs
１
６
８
０
０

１
２
８
０
１
０

６
４
０
１
１

CPUs+MICs
６
４ ３
２
１
２
１
１
０
１

Lecture Note

3. The Many-Core Accelerator and Examples of Its Use at IFERC-CSC

た．合計１２８ランクの１０２４×１０２４×７６８において３２ノードでN. Nakajima et al.

２．マルチ MIC によるネイティブモードの計算では１
２０ス

は MIC のみで実行することができなかった．また，
１６ノーレッドを使った計算が最速となった．

ドでは CPU+MIC で実行することができなかった．これら

３．MIC は CPU+ MIC の組み合わせであるヘテロジニアスの原因はMICのメモリ容量不足であり，MIC内に２ランク

計算において MPI でランクを作ることができるの

たてることができなかった．そのため，分割数を１２８からで，演算コアとして使用することができる．つまり，

２５６に変えた．

コア内でランクを振ることができるため，全体の並列

図１７に１０２４×１０２４×７６８において合計ランク数を２５６に固化効率を向上させることができる．

定した時，ホスト CPU，MIC，CPU+MIC を様々な組み合

４．ホスト CPU，ネイティブモード，およびこれらのハイ

わせで稼働させた時の速度比較を示す．表４は図１７におけブリッドの速度を比較すると，ホスト CPU が最速とる１つのCPU，１つのMICあたりのランク数の組み合わせなった．そして次は，シンメトリックモード，MICを示している．

ネイティブモードであった．１
０２４×１０２４×７６８におい

ランクの数の合計はどのパターンにおいても２５６となる．

てノードの数が３２，６４の時シンメトリックモードとホ全てのケースにおいて MIC 内のスレッド数は１２０に固定しスト CPU が同じ数のランクを使用する場合はシンメ

た．３２ノードと６４ノードでは CPU+MIC がホスト CPUの速トリックモードのほうがホスト CPU より速くなる．

度を超えた．しかし，１２８ノードを使った時はホスト CPU

ヘテロジニアススーパーコンピュータでシンメトリック

が最速となった．

モードが使える環境であるなら，上記のように少ないノー

３．
３．
４ 結論

ド数で高性能な計算を行うことができるということが明ら

速度比較の評価はホストCPUを使ったものや，MICのネかになった．

イティブモード，シンメトリックモード（ヘテロジニアス
計算）といった複数の計算方法で行った．Helios 計算機を参考文献

使用した場合の計算時間は以下の様な結論が得られた．

［１］中島徳嘉 他：プラズマ・核融合学会誌 91, 711 (2015).
［２］http://www.top500.org/lists/2015/11/.
［３］ジム・ジェファース，ジェームス・レインダース（訳：
すがわらきよふみ・エクセルソフト株式会社），インテ
ル!Xeon PhiTM コプロセッサ ハイパフォーマンス・
プ ロ グ ラ ミ ン グ（株 式 会 社 カ ッ ト シ ス テ ム，日
本，２
０
１
４）
．
［４］C.K. Birdsall and A.B. Langdon, Plasma Physics via Computer Simulation (Institule of Physics Publishing, Bristol
and Philadelphia, 1991).
［５］内藤裕志：プラズマ・核融合学会誌 74, 470 (1998).
［６］内藤裕志，佐竹真介：プラズマ・核融合学会誌 89, 245(2013).
［７］V.K. Decyk and T.V. Singh, Comput. Phys. Commun. 182,
641 (2011).
［８］V.K. Decyk and T.V. Singh, Comput. Phys. Commun. 185,
708 (2014).
［９］日本原子力研究開発機構 IFERC-CSC 利用委員会他編
集：プラズマ・核融合学会誌 92, 157 (2016).
［１０］佐竹信一 他： マルチ MIC＋マルチ CPU を使ったヘテ
ロジニアス計算を用いた MPI+OpenMP による円管内
乱流の直接数値計算 ，第２９回数値流体力学シンポジウ
ム講演論文集 (2015)12.15-17，九州大学，B08-1.
［１１］S. Satake et al., Fusion Sci. Tech. 68, 640 (2015). dx.doi.org/10.13182/FST14-956［１２］K. Hosaka et al., "DNS of a turbulent flow by MPI+
OpenMP for heterogeneous computing using multi-CPU
and multi-MIC", Proc. 28th International Conference on
Parallel Computational Fluid Dynamics, Parallel CFD
2016 (2016).
［１３］S. Satake et al., Phys. Fluids 18, 125106 (2006).

１．コードが MPI+OpenMP で書かれているとき，ホスト
CPU で並列計算のスレッド数を増加させることで，メッシュの数に比例する速度向上を得ることができ
る．

図１
７ ２
５
６ランクの場合のノード数変化による計算時間（１
０
２
４×
１
０
２
４×７
６
８）
．
表４

1024×1024×768（256 ranks）における CPU と MIC の組み
合わせ．

Nodes
data/CPUs
data/MIC０
data/MIC１

CPUs
MICs
CPUs+MICs
１
２
８ ６
４ ３
２ １
２
８ ６
４ １
２
８ ６
４ ３
２
２
４
８
０
０
１
２
４
０
０
０
１
２
１
１
２
０
０
０
１
２
０
１
２

６２１

Journal of Plasma and Fusion Research Vol.92, No.8 August 2016

核融合科学研究所 ヘリカル研究部 六ヶ所
研究センター 教授．今年５月から六ヶ所
赴任．主として MHD・新古典理論の研究に従事してきたが，IFERC 事業長兼任もほぼ６年となり，六ヶ所生活共々，得難い経験と人間関係を得た．事業の２０
１
９年末までの延長も決まり，新たな展開に期待している．

!!!!!!!!!!!!!!!!!!!!!!!!!!
とう

ひろ

し

山口大学大学院創成科学研究科教授．学部
担当は知能情報工学科．研究分野はプラズ
マのコンピュータシミュレーション．趣味
はハングルの勉強．本年度末で無事定年退職の予定．定年後は明るい下流老人として最高の人生を送りたいと思っています．

みや

と

なお

あき

宮 戸 直 亮量子科学技術研究開発機構・六ヶ所核融合研究所 上席研究員・BA 計画調整グルー
プサブリーダー．日本原子力研究所，日本
原子力研究開発機構を経て２０
１
６年４月から

現職．と書くと研究機関を渡り歩いているように見えるが，実質の所属は変わらず，所属機関の名前が十数年で２度変わっている．１回ならともかく，２回も変わるとは……こう
なると３回目もあるかも．専門は磁場閉じ込めプラズマの理論・シミュレーションで，最近はジャイロ運動論モデルの拡張などプラズマの基礎モデルの研究が中心．

!!!!!!!!!!!!!!!!!!!!!!!!!!
さ

たけ

しん

いち

佐 竹 信 一東京理科大学基礎工学部電子応用工学科
教授．専門分野／関心分野：数値熱流体力学／特に核融合炉冷却に関する MHD 乱流
熱伝達．

!!!!!!!!!!

!!!!!!!!!!!!!

ない内 藤 裕 志!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!

中 島 徳 嘉!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!!!!!

よし!!!!!!!!!!!!!!!!!!!!!!!!!!

のり!!!!!!!!!!!!!!!!!!!!!!!!!!

じま!!!!!!!!!!!!!!!!!!

なか!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!
き東京理科大学大学院基礎工学研究科電子応用工学専攻．専門分野／関心分野：シミュ
レーション工学／特にヘテロジニアス数値
計算．

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!

!!!!!!!!!!

かず

く

ぬぎ

とも

あき

功 刀 資 彰京都大学 物理工学系（大学院工学研究科）教授．専門分野／関心分野：伝熱工学，
特に沸騰現象の実験と解析／混相流工学／
数値熱流体力学／核融合炉工学（ブランケット，ダイバータ，安全性など）．

!!!!!!!!!!!!!!!!!!!!!!!!!!

さか!!!!!!!!!!!!

ほ保 坂 和 樹!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

６２２!!!!!!!!!!!!!!!!!!!!!!!!!!