J. Plasma Fusion Res. Vol.93, No.6 (2017)２７１‐２７７

研究・技術ノート

簡約化 MHD シミュレーションコードの自動チューニング
Automatic Tuning of Reduced MHD Simulation Code宮 川 大 和，TRETLER Rudolf，龍 野 智 哉
MIYAGAWA Yamato, TRETLER Rudolf and TATSUNO Tomoya
電気通信大学（原稿受付：２０１
６年１２月６日／原稿受理：２０１
７年４月１１日）

擬スペクトル法による２次元簡約化電磁流体コード rmhdper を，スレッド並列化により自動的にチューニングするシステムを開発した．本システムは並列スレッド数を自動で制御するものであり，ユーザーはグリッド数やマシンを変更するたびに最適スレッド数を選定する作業から解放される．本システムはコード中で変化するメモリアクセスパターンに応じて最適なスレッド数を個別に自動選択するため，コード全体で固定されたスレッド数を用いる場合と比べても実行時間が短縮されることがある．Helios で４０９６×４０９６グリッドを用いた場合，自動チューニングシステムは処理の種類によって１６，８と変化するスレッド数を自動的に選んでシミュレーションを行い，これはどの固定スレッド数を選んだ場合よりも高速であった．
Keywords:
magnetohydrodynamics, simulation, high performance computing, OpenMP, automatic tuning

１．はじめにの並列化に着目し，マルチスレッド並列処理を効率的に行自動チューニングは，計算機の性能向上が継続してきたう問題を取り扱う．

末に生まれた比較的若い研究分野である［１‐４］．その目的

マシンに搭載されている CPU コアを全て利用してマルにはアプリケーションを高速に動かすこと，処理のためのチスレッド並列処理を行った場合，メモリアクセスが律速消費電力を抑えることなど幾つかの視点があり，そのためとなる偏微分方程式ソルバーでは，メモリアクセスの衝突に配列の構成，アルゴリズムや利用ライブラリの選択，並が頻発して計算性能が低下することがある．したがって，列度の調整など，様々なレベルでの自動チューニングが考計算処理の種類に応じてあえてスレッド数を CPU コア数えられる．例えば Active Harmony［１］は独立した監視サーよりも少なくすることで，性能が向上する可能性がある．

バを立ち上げ，計算機のリソース状況に応じてライブラリ

最適なスレッド数は計算機の特性や処理内容，グリッドサ選択やパラメータ調整を動的に行うことで最適な高速化をイズなど複数の要因に依存するため，先見的に定めることめざす．また Autopilot［２］など，ネットワークを介して広は困難である．したがって，現実には試行により定めるこ

域に分散したシステムにも対応し，センサーからのデータとが多いわけであるが，この試行をプログラミングしておをファジィ推論に基づく決定プロセスを経て動作させるとき，自動的に定めることができればユーザーの負担を大幅いう壮大な研究もある．ここでは，並列度を自動的に調整に軽減することができる．本研究では，このような試行をすることで，数値シミュレーションコードを高速動作させプログラムにより自動的に行うものを考える．

ることを考える．

本研究では，コードの可読性と確実性を重視して設計さ近年の計算機では CPU の動作周波数は頭打ちしており，れた，単一ノードで動作する簡約化 MHD コード rmhdper

代わりに複数のコアを搭載することによる性能向上が図ら［５］の高速化を考え，計算の種類に応じて最適なスレッドれている．大型計算ではそのような複数のCPUコアを搭載数を自動的に選択する自動チューニングシステム（Auto-した計算機を複数台（ノード）用いることにより，ノードmatic Tuning System: ATS）を開発した．このシステムに内，ノード間の階層的な並列化が必要となってきた．これ

より，ユーザーは最適なスレッド数を試行により決定する

らノード内，ノード間の異なる階層における並列処理でプロセスから解放される．さらに，複数のマルチスレッドは，共有メモリ環境で用いられるマルチスレッド並列処理並列処理部に異なるスレッド数を用い，最適な並列計算をと，分散メモリ環境にも対応できるマルチプロセス並列処自動で行うことができるようになった．この手法は，分散理の二種類を効率的に用いる必要がある．近年はプロセスメモリと共有メモリを併用したさらに大型のハイブリッドとスレッドの使い分け方も，シミュレーションを行う研究

並列計算コードにおいても，共有メモリ並列部分にそのま者を悩ませる新たな負担となっている．ここではノード内ま用いることができる．

University of Electro-Communications, Chofu, TOKYO 182-8585, Japancorresponding author’s e-mail: tomo.tatsuno@uec.ac.jp

２７１!2017 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.93, No.6 June 2017

まず次節で本研究で用いる数値コード rmhdper の概要"!
#*)%について述べたあと，第３節で OpenMP による並列化の実装方法を示す．第４節で我々の開発した ATS について述# '#"''"'!
% '&%''&'（３）を表す．プラズマの流速ベクトル "'%!
&(と磁場ベクトル%を '方向の単位 ベ ク ト ル と し て そ れ ぞ れ!'%!
&(は，#べ，第５節でまとめと今後の展望について議論する．
本研究では，電気通信大学の研究室にあるデスクトップ%%と表される．
，!%&&##"%&###機（ローカルマシンと呼ぶ）と，国際核融合エネルギー研
究センターの計 算 機シミュレーションセンター（IFERC-CSC）

本研究では，周期境界条件の下で方程式（１）−（２）の初にある大型計算機 Helios を利用した［６］．ローカルマシンの期値問題を解くオープンソースのコード rmhdper［５］を使CPU はクアッドコアの Intel Core i7 920，Nehalem 2.67 GHz用する．空間の離散化には擬スペクトル法［８］を用いるが，であり，OS は Linux Mint 17.1 Rebecca，コンパイラは In-

そ の 際 の 高 速 フ ー リ エ 変 換（Fast Fourier Transform:

tel Fortran Compiler 14.0.1 を用いた．Helios の CPU は１６コFFT）は FFTW ライブラリ［９，１０］を利用する．時間の離散

アの Intel Xeon E5-2680，Sandy-Bridge 2.7 GHz であり，OS化には線形項について２次の陰的後退差分法［１１］を，非線は Red Hat Enterprise Linux 6.4 Santiago，コンパイラは In-形項について３次の Adams-Bashforth 法（AB3）［１２］を使tel Fortran Compiler 15.0.2 を 用 い た．ま た コ ー ド の プロ用する．ただし，時間ステップの１ステップ目には線形項ファイル取得には，GNU プロファイラ gprof を用いた．

および非線形項ともに Euler 法を使用し，２ステップ目には非線形項に２次の Adams-Bashforth 法を使用する．
２．
２ 変更前コードのプロファイル

２．簡約化 MHD ソルバー rmhdper の概要
２．
１ 支配方程式と数値解法

変更前の rmhdper では，処理性能にはあまり注意を払わ

MHD 方程式に非圧縮性と２次元性を仮定して簡約化れていなかった．変更を加える前に，gprof によるプロファ

し，Alfvén 時間で規格化された２次元の簡約化 MHD 方程イル結果（一部抜粋）を表１に示す．また比較のため，以式［７］降に述べる変更をすべて行った最終的なプロファイル結果も同時に表２に示した．計測はローカルマシンで行い，グ'"#")
#!
"#*
&!
"&*
""#，%)
"%
'
$

（１）

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

（２）

リッド数を '
，ステップ数を５００とし!%!
!
!&("!#$"!#$(%'た．gprof ではサブルーチンごとに呼び出し回数，処理時間などが計測されるが，表１の各列はそれぞれ以下を表す：を考える．ここで#'
%!
&(と&'
%!
&(はそれぞれ流れ関数と磁% time

束関数を表し，定数 %と $はそれぞれ動粘性係数と電気抵全体の総実行時間に対する，サブルーチンの総実行時間の占める割合，抗率を表す．また "%&"は２次元のラプラシアン，)!
$$*

cumulative seconds プロファイルリストにおいて，該当は Poisson 括弧列を含め上位に示されている全てのサブルーチンの累

表１００．
変更前の rmhdper コードのプロファイルの一部．グリッド数(Nx, Ny) ＝ (2048, 2048)，ステップ数５%
time
37.22
15.18
9.84
9.01
5.96
4.66

表２

self
seconds
165.47
67.50
43.73
40.06
26.52
20.70

calls
5068

self
s/call0.03

total
s/call
0.03

1500
500
1002

0.03
0.05
0.02

0.10
0.42
0.02

name
fft_mp_ktox__intel_sse2_rep_memset_intel_ssse3_rep_memcpy
four_mp_poisson_bracket_
tint_mp_advect_
fft_mp_xtok_

変更後の rmhdper コードのプロファイルの一部．グリッド数(Nx, Ny) ＝ (2048, 2048)，ステップ数５
０
０．

%
time
32.85
11.79
10.32
10.30
6.23
1.12

cumulative
seconds
165.47232.97276.70316.76343.28363.98
後略

cumulative
seconds
70.91
96.37
118.65140.89154.34
中略
199.39後略

self
seconds
70.91
25.46
22.28
22.24
13.45

calls
4316

self
s/call0.02

total
s/call
0.02

1072
536
1076

0.02
0.04
0.01

0.06
0.20
0.01

2.42

name
fft_mp_ktox__intel_ssse3_rep_memcpy
four_mp_poisson_bracket_
tint_mp_advect_
fft_mp_xtok__intel_sse2_rep_memset

２７２Research and Technology Note

Automatic Tuning of Reduced MHD Simulation Code

Y. Miyagawa et al.

回バージョン３．
３．
４において使用したプラン作成関数は，計実行時間（秒），self seconds サブルーチン単独の総実行時間（秒），実数から複素数配列への dfftw_plan_dft_r2c_2d と複calls サブルーチンが呼び出された総回数，素数から実数配列への dfftw_plan_dft_c2r_2d であself s/call サブルーチン自身の一回の呼び出しにかかる

る．プラン作成の際１６バイト境界にアラインされた配列も指定し，変換の実行には，dfftw_execute を２次元の１平均実行時間（秒），total s/call 一回の呼び出しあたりの，自分自身とその内変換あたり１度だけ実行する．キャッシュ効率を考慮して部で呼び出す子サブルーチンの実行時間の和の平均１次元変換を組み合わせることで更なる高速化を計ること（秒），も考えられるが，今回は並列度の変更による自動チューニ

nameングのみを考え，アルゴリズムの変更は含めないこととしプロファイルを取得したサブルーチンの名前．

た．また，プラン作成にあたって使用した最適化オプショ
ンは FFTW_PATIENT である．

プロファイル結果は self seconds により降順に並べられ，表１ではその上位６つを示した．表中の空白は，コンパイさらに呼び出し回数削減による高速化を２つ行った．コーラの内部関数など，呼び出し回数が決定されなかったこと

ド作成時の確実性のために，実行時に割り当てた配列は一を示す．各サブルーチンの内容は以下のとおりである：

旦全要素にゼロ代入を行って初期化し，その後に必要な分だけ代入などの操作が行われていた．サブルーチン内で呼び出しごとに割り当てられる配列ではこの操作が多数回にfft_mp_ktox_ 波 数 空 間 に お け る 配 列 を 実 空 間 へ 逆及び，それ が 表１に お け る_intel_sse2_rep_memset

フーリエ変換（IFFT）
する，の肥大化につながっていると考えられる．実際の代入がす

fft_mp_xtok_ 実空間における配列を波数空間へフーべての要素に及ばないこともあるため，いくらかのゼロ代リエ変換（FFT）
する，
four_mp_poisson_bracket_

入は必要であるが，必要最低限に限ることで代入回数を削

FFT，IFFTや微分演算減した．

関数を呼び出し，Poisson 括弧を計算する，
tint_mp_advect_ Poisson括弧やAB３を計算する関数

呼 び 出 し 回 数 削 減 の２つ 目 と し て，IFFT の 回 数を呼び出し，非線形項の時間ステップを進める，を 削 減 し た．表１か ら ス テ ッ プ 数 が５００の 場 合，_intel_sse2_rep_memset，_intel_ssse3_rep_memcpy

four_mp_poisson_bracket_の変数への値の代入，コピーを行うIntelコンパイラの内回

数が １５００ 回，fft_mp_ktox_の回数が約５０００回であることがわかる（余剰分は結果の診断などに用いられる６８回）．つまり，表１で部ルーチン．

は１時間ステップあたりfour_mp_poisson_bracket_が
表１より，プロファイルが取得されたサブルーチンの中３回，fft_mp_ktox_が１０回呼び出されていたことにな

で fft _ mp _ ktox _， four _ mp _ poisson _ bracket _ ，る．本コードで使用する方程式（１），（２）では，Adams-tint_mp_advect_の順に総実行時間が大きいことがわか

Bashforth 法を用いる場合，１時間ステップあたり３つのる．高 速 化 を 図 る 場合，サ ブ ル ー チ ン の呼 び 出 し 回 数

Poisson 括弧を１度ずつ評価する必要がある．これらを（calls列）を減らしたり，一回の呼び出しあたりの平均実行

FFT !，IFFT !

時間（self s/call 列，total s/call 列）を短くすることで高速#!
"#('#!

!!

を陽に用いて表せば!!
$
$
"
#%#
#$#"#%!"&&! %!!
!!
$
$
"
#$#&
#%#"#!"&−! %! %
，化を実現できる．各サブルーチンの calls 列と total s/call列 に 注 目 す る と，fft_mp_ktox_は 呼 び 出 し 回 数!!

$
$
"
#%$
#$#"$
&&%!"
! %!!
!!
"$$
"
#$$&
#%# $&!"−! %! %，$!
"$('#!

（calls 列）が５
０６８回と最も多いために時間がかかっていること が わ か る．一 方，four_mp_poisson_bracket_や!!

!!
$$
"
#%#"
#$$&%&! %
!!
!!
$$
"
#$#&"
#%$
&，
−! %
! %

#!
$('#!

tint_mp_advect_は内部に留まる平均実行 時 間（totals/call 列）が長いために時間がかかっている．以下では，

（４）

!!

（５）

!!

（６）

$#
$#，$$&である．各
#&
##!%
#!%
物理量は基本的に波数空間でメモリに保存されており，表１の上位６つに関する領域について高速化を図ることととなる．ただしここでする．

Poisson 括弧の計算でのみ実空間に変換して，この後再び３．rmhdper の高速化

波数空間に変換して時間発展計算を行う．
（４）‐（６）では３．
１ 逐次実行レベルの高速化

１２回の IFFT（!

!!

）が現れるが，これをいくらか削減して前節で挙げた表１の上位６つのサブルーチンについて，１０回で行っていたのが表１である．つまり，Poisson 括弧は逐次実行のままで高速化できないか検討した．まず FFT

#!
!
'(#!"
$!であることを利用し，!を先に計算しておく．
four_mp_poisson_bracket_をベクトル !とスカラーに用いられるライブラリ FFTW のバージョンは２．
１．
５を利た．３．
３．
４は単独CPUの使用でも２．
１．
５よりいくらか高速化

量 !を受け取るように設計することで，（４）と（６）で重複!!
!!
$$"
#%#"
#$#&&の計算を１回ずつ削減していする! %，! %されているが，新たに OpenMP にも対応しており，次節でたわけである．

用していたが，これを３．
３．
４も利用できるように 変更しところが，（６）に出てくる残りの!の並列化に利用される．FFTW は，実際に FFT を行う前!!
$$"#$$"#%$&&%，! %!!

は既に（５）で計算したものである．これらも同時に再利用に初期化段階で変換のプランを作成する仕組みである．今２７３

Journal of Plasma and Fusion Research Vol.93, No.6 June 2017

することで，（６）は four_mp_poisson_bracket_を呼び最後に，tint_mp_advect_から呼び出される３次の出すことなく計算できることがわかる．したがって１時間

Adams-Bashforth 法の計算領域 ab3 では，ブロック分割とス テ ッ プ あ た り，four_mp_poisson_bracket_を２４］，およびそれらの中間と言えるブ
サイクリック分割［１回，fft_mp_ktox_を８回呼び出すだけでよい．

ロックサイクリック分割を試した．これら３手法で並列化

３．
２ OpenMP による並列計算の実装

効率に有意な違いは見られなかったため，本研究では最後に試したブロックサイクリック分割を用いている．

本節では，表１に挙げられたコード中の処理時間が大きい３つの領域について，個別にスレッド並列処理を実装す

４．自動チューニングシステム（ATS）

る．

４．
１ 背景

ま ず 表１の サ ブ ル ー チ ン の 中 で，fft_mp_ktox_やfft_mp_xtok_は内部でFFTWライブラリを呼んでおり，マルチスレッド並列処理による並列計算を行う場合，各FFTW3 のスレッド並列処理を利用することができる．こマシンの最大利用コア数を利用スレッド数として固定しての並列処理領域を fft と呼ぶことにする．FFTW3 の変換計算［図３（a）］することが多い．しかし，スレッド数が大きプラン作成では，プラン作成前にマルチスレッド処理命令い場合，メモリアクセス量が演算量を上回ったり，並列処文を追加することで並列化が可能となる．図１のように!$理のための分配，分割および結合といったオーバーヘッドから始まる OpenMP 指示文を利用して，FFTW3 のプランが大きくなる可能性がある．そのため，共有並列計算に利作成前にマルチスレッド処理命令文を追加した．fft の並用可能な搭載コア数が必ずしも最適なスレッド数とは限ら列処理に使用するスレッド数は nthreads_fft で指定さない．さらに，コード内にある様々な計算領域によって演れる．ただし，iret は算量やメモリアクセスのパターンが異なるため，最適なスFFTW3 のスレッド初期化時のエレッド数も異なることが考えられる．

ラーチェック変数であり，エラーが発生した場合は０が，そこで，スレッド数を実行時に自動的に調節する自動正常な場合は１が代入される．
次に，IFFT におけるエイリアシング［１３］に関連した配列コピー領域 ktox の並列化を考える．図２にサブルーチン ktox の 一 部 を 示 す が，こ れ は 波 数 空 間 に お け る
complex 型の２次元配列 fa を受け取って，real 型の２次元配列 fyx を返すものである．図２に示されたとおり，サブルーチン fft_mp_ktox_全体に使用するスレッド数をnthreads_ktox で指定して並列処理実行指示文を与えた．ループ構文と構造化されたブロックに対しては do 指
示文と workshare 構 文 に よ っ て 並 列 処 理 命 令 を行い，FFTの実行部分に対してはmaster構文を使用することでマスタースレッドのみが dfftw_execute関数をコールするようにした．つまり FFT のスレッド並列化は，スレッドの立ち上げを含め全てがFFTW側で処理されている．またmaster 構文はスレッド間の同期を自動で行わないため，IFFT 実行後に barrier 指示文によって明示的に同期を行っている．

（以下，プランの作成）
図２
図１

fft に関するスレッド並列計算．

２７４

ktoxの並列領域の実装箇所．最初の!$OMP PARALLELの行は，紙面の都合上折り返されている．

Research and Technology Note

Automatic Tuning of Reduced MHD Simulation CodeY. Miyagawa et al.

チューニングシステム（ATS）を開発した．このシステム

る順序は，計算時間の大きい順に fft，ab3，ktox とした．

は，初期化時にスレッド数を変化させながら数回のフルタ

また，比較したいスレッド数の最大値は２の乗数のうちかイムステップを計算し，経過時間の計測結果をコード内でらインプットファイル内で設定することができ，未設定の比較することで最適なスレッド数を自動で選択しようとす場合は環境変数 OMP_NUM_THREADS の値が用いられる．自るものである．このシステムによりコード中の各領域にお動チューニングに利用する時間ステップ数は，各グリッドけるスレッド数が個別に最適化［図３（b）］され，計算中の数で１回の計測が１秒程度に収まるよう，グリッド数に応スレッド数変更も含めた最適な並列計算を自動で行うこと

じて２から１００の間から適当に選んだ．
図４は 実 際 に rmhdper コ ー ド を 実 行 し た 時 の 自 動が期待できる．
４．
２ システムの概要

チューニング結果の出力例である．ただし，実行はローカ

前節で決定した３箇所の並列領域 fft，ktox，ab3 につルマシンで行い，グリッド数は "!
!$#!#!
#!"#!"#を用!"いて自動チューニングを行った．本コードでの ATS の処

い た．図４で は，各 並 列 処 理 領 域 ご と に ス レ ッ ド 数

理手順は以下のとおりである：

threads と計測時間 time が並んで出力されている．時間計測の結果，各並列処理領域でのスレッド数は fft が８ス

１．イ ン プ ッ ト フ ァ イ ル 中 で 設 定 さ れ た 論 理 変 数

レッド，ab3 が２スレッド，ktox が２スレッドに決定されautotune の値が真の場合，コードの初期化時に自動ていることがわかる．

チューニングが開始される，

４．
３ 性能評価

２．まず ab3，ktox のスレッド数を１に固定し，fft のス

簡約化 MHD コード rmhdper について，レッド数のみを１，２，４，…と変化させながら実行時間を計測し，比較する．最も実行時間が短い時のスレッド

１）初期状態のコード（Before），数を fft の本計算でのスレッド数とする，２）固定スレッド数を用いた並列化（"threads），
３）ATS を利用した場合，３．決定した fft のスレッド数はそのまま，ktoxのスレッド数を１に固定し，ab3 のスレッド数を１，２，４，…と変における計算時間の比較を行い，性能を検証した．

化させながら実行時間の計測，比較を行う．最も実行計測には研究室のローカルマシンとスーパーコンピュー時間の短いスレッド数を ab3の本計算でのスレッド数

タ Helios を利用した．計測結果はすべて１ステップあたりとする，の時間［秒/１ステップ］に焼き直すが，本計算のステップ４．決定した fft と ab3 のスレッド数はそのまま固定し，数はグリッド数ごとに異なる．グリッド数が６４×６４の時はktox のスレッド数を１，２，４，…と変化させながら実行時間の計測，比較を行う．最も実行時間の短いスレッド数を ktox の本計算でのスレッド数とする，５．各領域で使用するスレッド数が決定され，本計算を行う．
上記の計算手順において，並列領域のスレッド数を決定す

図３

固定スレッド並列時と ATS 利用時との並列処理領域のスレッド数の違い．（a）
固定スレッドによる計算実行時，各並列領域におけるスレッド数は同一となる．（b）
ATS を利用した計算実行時，各並列領域におけるスレッド数は異なる最適化がされる．

（以降，本計算アウトプット）
図４２７５

ローカルマシンにおけるグリッド数(Nx, Ny) ＝ (512, 512)の自動チューニングの様子．

Journal of Plasma and Fusion Research Vol.93, No.6 June 2017

ステップ数を８００００回，１２８×１２８の時は２００００回，２５６×２５６を上回らないように並列数を自動的に調整している．の 時 は ５０００回，５１２×５１２ の 時 は ２０００回，１０２４×１０２４，４．
４ 変更後コードのプロファイル

２０４８×２０４８および４０９６×４０９６の時は５００回に調整した．シ

変更後の rmhdper コードに対して，gprof を使って変更

ミュレーションデータの出力は，計測時間には含めないこ

前と同条件で再度プロファイルを取得した結果の一部をととした．性能評価を行うにあたり，実行時間の決定は，表２に示す．表１と比べると，関数 fft_mp_ktox_，four

同条件の計算を５回ずつ行い，それぞれの計算時間を計測_mp_poisson_bracket_，tint_mp_advect_すべてに後，その中の最高値と最低値を除いた残り３つの値の平均おいて，１回の呼び出しに対する自分自身と子サブルーチ値を用いた．

ンの平均時間（total s/call）が短くなっていることがわか

図５に計算時間の比較結果を示す．図５（a）はローカルる．これは，マルチスレッド並列処理を施し上記３つの関マシンでグリッド数１０２４×１０２４のとき，図５（b）は Helios数に関わる領域を高速化したためだと考えられる．また，でグリッド数４０９６×４０９６のときの１時間ステップ当たりにfour_mp_poisson_bracket_の呼び出し回数（calls）が掛かる時間を示す．変更前（ラベル B）と固定１スレッド

１５００から１０７２に減少したのは第３節で述べたように，Pois-（ラベル１）はいずれも逐次実行であるが，ケース（b）におson 括弧の演算回数を１ステップあたり２回呼び出しを行いて 11.2秒から4.44秒へ短縮されており，これは３．
１節で述うように変更し，また自動チューニングの初期計測によっべた高速化による．自動チューニング（ラベル A）が最もて７２回 呼 び 出 さ れ た た め で あ る．こ れ に 伴 い，計算時間が短く［ケース（b）では 1.72 秒］，どの固定スレッドよりも高速化されていることがわかる［ただしケース（a）では４スレッド（ラベル４）とほぼ同等］．
さらに比較するグリッド数の範囲を広げた結果を図６および７に示す．図６はローカルマシン で の 実 行 結果を，図７は Helios での実行結果を表す．図６，７において ATS利用時の計算時間は "印付き実線で示されており，常に他の固定スレッド利用時より短い，もしくは同等となっていることがわかる．一方，固定スレッド利用時は，例えばHelios の１６スレッド利用時（!印付き破線）を見ると，グリッド数が大きいときは計算時間が短いが，グリッド数が小さくなると計算時間が比較的長いことがわかる．
このように，搭載CPUコア数を使用スレッド数として固図６定した場合の計算時間が最も短いわけではない．ATS が実

ローカルマシンでの計算時間比較図．

際に選んだ最適なスレッド数のうち，特徴的なものを表３に示す．実行のタイミングによって最適スレッド数は多少変化することがあるが，ローカルマシン，Helios の双方でグリッド数が大きいほどスレッド数も大きくなる傾向を示
す．３つの並列化部分のうちでは fft が最もスレッド数が多く，ab3 や ktox は同程度で，fft に比べて少なめのスレッド数が使われている．これは，並列計算のためのオーバーヘッドが発生しているためであり，ATS を利用した場合はオーバーヘッドが並列化によって得られる高速化効率

図７

ATS が選んだグリッド数別の特徴的なスレッド数．

local

表３計 算 時 間 比 較 図．（a）
１０２４×１０２４ロ ー カ ル マ シ ン，（b）
０９６スーパーコンピュータ Helios．横軸のラベルは，
４０９６×４
B は変更前の逐次実行コード（Before），数字は固定されたスレッド数，A は ATS を表す．

Helios

図５２７６Helios での計算時間比較図．

grid
１２８2５１２2２０４８2２５６2１０２４2４０９６2

fft
１８８１６１６１６ab3１２２１８８ktox２２２２８８Research and Technology Note

Automatic Tuning of Reduced MHD Simulation CodeY. Miyagawa et al.

fft_mp_ktox_の呼び出し回数（calls）も減っている．さら

介氏には，システムの開発にあたり多くの有益な助言を得に，省略可能な配列への代入文を削除したことにより，た．名古屋大学の前山伸也氏には文献［１４］を紹介いただい_intel_sse2_rep_memset の実行時間がかなり減少した．また匿名の査読者には，本原稿の改善につながる多数た．の有益なコメントをいただいた．

５．まとめ

参考文献

本研究では，擬スペクトル法による簡約化 MHD コード［１］C. Tapus et al., the IEEE/ACM SC2002 Conf. (2002).
［２］R.L. Ribler et al., Future Gener. Comput. Syst. 18, 175(2001).
［３］弓場敏嗣：電気通信大学紀要 20, 1 (2007).
［４］須田礼二：情報処理学会研究報告 2011-HPC-129 (14)(2011).
［５］The Gyrokinetics Project, https://sourceforge.net/p/
gyrokinetics/code/HEAD/tree/rmhdper/［６］日本原子力開発機構 IFERC-CSC 利用委員会，IFERC
事業チーム，IFERC-CSC ホームチーム：プラズマ・核融合学会誌 92, 157 (2016).
［７］H.R. Strauss, Phys. Fluids 19, 134 (1976).
［８］J.P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed.
(Dover, Mineola, 2001).
［９］M. Frigo and S.G. Johnson, Proc. IEEE 93, 216 (2005).
［１０］FFTW Homepage, http://www.fftw.org/［１１］C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hall, Englewood Cliffs, 1971).
［１２］J.D. Hoffman, Numerical Methods for Engineers and Scientists, 2nd ed. (Marcel Dekker, New York, 2001).
［１３］S.A. Orszag, J. Atmos. Sci. 28, 1074 (1971).
［１４］渡邉智彦，井戸村泰宏：プラズマ・核融合学会誌 89,171 (2013).
［１５］R. Numata et al., J. Comput. Phys. 229, 9347 (2010); ibid.
245, 493 (2013).

rmhdper を簡易に高速化するため，スレッド並列度を自動で調節するチューニングシステムを開発した．このシステ
ムにより，複数のマルチスレッド並列処理領域におけるスレッド数が個別に最適化され，自動で最適な並列計算を行うことができるようになった．利用者はグリッドサイズや計算機を変更するたびに最適なスレッド数を試行により選定する作業から解放され，さらにどの固定スレッド数よりも高速な実行が実現される．
グリッド数 "!
!"!
!##"!$#"!$##のシミュレーション!"を Helios で行ったところ，１ステップあたりの計算時間は，固定１スレッドの 4.44 秒から，８スレッドで 2.08 秒，１６スレッドで 2.43 秒，自動並列処理（１６，８スレッド混在）で 1.72 秒へと短縮された．
今回作成した自動並列システムは MPI との併用が可能である．AstroGK［１５］などの MPI 並列コードに今回の手法を実装し，ハイブリッド並列にも応用することができるだ
ろう．

謝

辞

本研究の大部分は，国際核融合エネルギーセンターの計算機シミュレーションセンター（IFERC-CSC）における大型計算機 Helios を用いて行われた．兵庫県立大学の沼田龍２７７