渡辺宙志
千葉修一
※ 本稿はメニーコア時代のアプリ性能検討WG 成果報告書にて公開された内容を整形したものです。
HPC分野においては、FortranないしC/C++言語が広く使われている。Fortran、C言語、C++言語のいずれも活発に仕様が改定され、現在も発展し続けている言語ではあるが、数値計算向けの「普段使いの言語」としては、おそらくPythonやJuliaなどの方が広く使われていると思われる。特にPythonは、昨今の機械学習ブームの牽引役として広く使われているようである。今後、PythonやJuliaといった言語によるスパコン利用も増えていくと思われるし、それにとどまらずRustのような次世代言語も使われるようになって欲しいと期待もするのだが、スパコンを使うためにFortran/C/C++が必要となる現状は当面続くであろう。
さて、スパコンを使うからには、プログラムは高速に実行されて欲しい。そのためには優秀な最適化機能を備えたコンパイラがあることが望ましい。ところが、FortranやC言語はともかく、言語仕様が巨大でプログラムの記述自由度が高いC++言語はコンパイラから見ると最適化が困難な言語である。スパコンユーザは、普段のコード開発をx86系のPC上で、GCC/clangもしくはインテルコンパイラを用いて行っていることが多い。x86上におけるGCCやclangの最適化機能は高度であり、かなり抽象的なレベルで最適化を行う能力を持っている。また、インテルコンパイラはSIMDベクトル化や、大域最適化(Link Time Optimization, LTO)に優れている。スパコンユーザはそういった「優秀」なコンパイラを使い慣れているため、富士通C++コンパイラの挙動に戸惑う場合がある。本WGでは、渡辺が2017年7月25日に富士通のC++コンパイラの問題点と改善の要望を挙げた。また、関連して2018年1月24日にマルチスレッド環境下におけるメモリ管理の問題について報告した。コンパイラについての問題については、富士通の千葉が2017年10月20日と2018年10月19日の二回にわたり、富士通のコンパイラ改善の現状と今後の見通しについて発表した。メモリ管理の問題については、2019年2月25日に富士通より報告があったが、こちらの内容については別の報告書にまとめ、本稿では主にコンパイラに対する要望と改善内容について述べる。なお、本WGにおいて渡辺は分子動力学法における力の計算のAVX2、AVX-512を用いたSIMDベクトル化手法とその性能向上についても報告を行ったが、本稿において詳細には触れない[1]。ポスト京アーキテクチャはSIMD幅が512ビットになることが発表されており、今後ますます性能におけるベクトル演算の重要性は増していくと思われる。コンパイラによる自動ベクトル化も重要なトピックと思われるが、本稿ではC++言語に特有の最適化にフォーカスする。
C++言語のようなクラスベースのオブジェクト指向言語においては、いくつかの変数や関数をまとめてクラスとして管理する。 渡辺が開発した分子動力学法コードMDACP[2]では、座標と運動量をまとめて
class Variables{
public:
double C0;
double q[N][D],p[N][D];
};という変数クラスを作って管理していた。Nが原子数、Dが次元を表し、座標q及び運動量pをそれぞれ多次元配列として宣言してある。また、C0はカットオフに使われる定数である。このVariableクラスのインスタンスを受け取り、力の計算をするコードを書いていたのだが、力の計算だけを抜き出したコードは最適化され、満足いく性能が出るにもかかわらず、それとほぼ等価と思われる本番コードは最適化されないことに悩んでいた。本件について調査した結果、コンパイラの最適化機能がクラス内のメンバ変数の定義順序に依存することがわかった。例えば、以下のようなトイコードを考えてみる。本コードでは変数C0は使われない。
void
calcforce(Variables *vars){
double (*q)[D] = vars->q;
double (*p)[D] = vars->p;
for(int i=0;i<N-1;i++){
const double qx = q[i][X];
const double qy = q[i][Y];
const double qz = q[i][Z];
double px = p[i][X];
double py = p[i][Y];
double pz = p[i][Z];
for(int j=i+1;j<N;j++){
double dx = q[j][X] - qx;
double dy = q[j][Y] - qy;
double dz = q[j][Z] - qz;
double df = (dx*dx + dy*dy + dz*dz);
px += df*dx;
py += df*dy;
pz += df*dz;
p[j][X] -= df*dx;
p[j][Y] -= df*dy;
p[j][Z] -= df*dz;
}
p[i][X] = px;
p[i][Y] = py;
p[i][Z] = pz;
}
}これはVariableクラスのインスタンスを受け取り、距離から定まる適当な力を計算して運動量に足し込むコードである。
冒頭で座標や運動量のポインタを定義して代入しているが、これは(少なくとも京稼働当初の富士通C++コンパイラでは)こうしないと最適化がかからなかったためである。
さて、このコードをそのままコンパイラに食わせるとソフトウェアパイプライニング(Software pipelining, SWP)がなされない。
「京」や「FX10」は、豊富なレジスタを活かして力任せにループをソフトウェアパイプライニングすることが前提のアーキテクチャとなっており、ソフトウェアパイプライニングが効かないと性能が出ないことが多い。しかし、ここでVariableクラスの宣言順序を入れ替えてみる。
class Variables{
public:
double q[N][D],p[N][D];
double C0; // 宣言を後ろにずらした
};すると先程の力計算のループにソフトウェアパイプライニングが適用され、実行性能が倍以上向上した。
今回のケースをユーザから見ると、「使われないクラスのメンバ変数の定義位置により、プログラムの性能が倍以上変わる」という奇妙な結果となる。
後に、restrict修飾されたポインタ間の最適化処理に問題があったのが原因と報告された。
このようにクラスとして保持されたデータのポインタを受け取る、という処理はC++のコンテナでは頻出するイディオムであるため、
同様な理由で性能が出ないユーザコードが他にもあることが想定される。この件については今後改善を検討する旨が報告された。
C++言語においては、生の配列を使わずにstd::vectorといったコンテナを使うのが一般的である。
コンテナに対するループ処理も、ループカウンタではなくイテレータを用いる。コンパイラとしてGCCやインテルコンパイラを使っている場合、コンテナを使っても生の配列を使った場合に比べて特に性能は劣化しない。しかし、富士通C++コンパイラはSTLコンテナを使うとループの最適化がうまくできない場合がある。
先程と全く同様なコードを、std::vectorを使って書き直してみよう。
const int N = 20000;
struct pos{
double x,y,z;
};
void
calcforce(std::vector<pos> &q, std::vector<pos> &p){
for(int i=0;i<N-1;i++){
const double qx = q[i].x;
const double qy = q[i].y;
const double qz = q[i].z;
double px = p[i].x;
double py = p[i].y;
double pz = p[i].z;
for(int j=i+1;j<N;j++){
double dx = q[j].x - qx;
double dy = q[j].y - qy;
double dz = q[j].z - qz;
double df = (dx*dx + dy*dy + dz*dz);
px += df*dx;
py += df*dy;
pz += df*dz;
p[j].x -= df*dx;
p[j].y -= df*dy;
p[j].z -= df*dz;
}
p[i].x = px;
p[i].y = py;
p[i].z = pz;
}このコードはソフトウェアパイプライニングされず、性能が出ない。
可変長配列であることが悪いのかと、固定長配列コンテナであるstd::arrayを使っても最適化されず、
引数に生配列を渡すと最適化が効き、妥当な性能が出る。なお、GCC、インテルコンパイラともにstd::vector、std::arrayを使っても生配列を使った場合と遜色ない性能がでることは付記しておきたい。
後に、STLの内部で使われているvolatile変数が最適化を阻害していることが報告された。こちらについても改善を検討中とのことである。
std::vectorその他、要素を多数保持するコンテナクラスについて、その要素全てについてなんらかの処理をする、といったプログラムは頻出する。例えば、std::vector<int> vが保持する要素全てを2倍したい、という時、従来のfor文を使ったプログラムでは以下のように書けるであろう。
for (unsigned int i = 0; i < v.size(); i++) {
v[i] *= 2;
}これは、C++11から導入されたranged-base forを用いて以下のように書ける。
for (auto &t : v) {
t *= 2;
}さて、通常のforを使ったコードとranged-base forを使ったコードは等価なのであるから、同じように最適化処理がなされて欲しい。
以下のようなコードを書いてみよう。
const int N = 10000000;
double a[N];
void
myabs(void){
for(auto &v: a){
if(v < 0.0)v = -v;
}
}単純に、グローバル配列aについて、その絶対値を取る関数である。このコードをコンパイルさせると、
SIMD化もソフトウェアパイプライニングも適用されない。通常のforを使ったコードに書き換えてみる。
const int N = 10000000;
double a[N];
void
myabs(void){
for(int i=0;i<N;i++){
if(a[i] < 0.0)a[i] = -a[i];
}
}こちらはSIMD化もソフトウェアパイプライニングも適用されており、アセンブリを見てもfneg,sやfcmplted,sといったSIMD命令が発行されていることが確認できる。
後にこれはループの解釈、とくに終了判定がカウント型(不等式による判定)であるか、終端判定型(等式による判定)であるかの違いが最適化に影響していること、後に改善を検討することが報告された。
C++では、データのまとまりを構造体として扱うことが多い。 また、演算子オーバーロードを用いて、構造体やクラスをあたかも基本型として演算したりするなど、柔軟な記述ができる。 以下は、二つの整数の組を要素にもつベクトルを表現することを意図した構造体である。
struct IntPair{
int i,j;
};
IntPair add(IntPair x, IntPair y){
IntPair z = {x.i + y.i, x.j + y.j};
return z;
}
void
func(int &x, int &y){
IntPair a = {1,2}, b = {3,4};
IntPair c = add(a,b);
x = c.i;
y = c.j;
}IntPairは整数を二つまとめた構造体で、IntPair addはその演算子オーバーロードを模した関数である。
関数func内では、二つのペアの和を計算し、その要素を参照の形で返却している。
さて、この関数の実行結果は常にx = 4, y = 6となるが、最近のコンパイラはそれを認識して、
無駄な構造体を作らず、値を即値で返してしまう。例えば以下はg++ -O3 -Sの実行結果である。
func(int&, int&):
movl $4, (%rdi)
movl $6, (%rsi)
ret4と6を即値で返していることがわかる。しかし、富士通コンパイラは関数funcについて以下のようなアセンブリを吐く。
func(int&, int&):
add %sp,-240,%sp
mov 1,%xg3
mov 2,%xg4
stw %xg3,[%sp+2255]
mov 3,%xg5
mov 4,%xg6
mov 3,%xg7
mov 4,%xg8
mov 1,%g1
mov 2,%g3
add %g1,3,%g2
add %g3,4,%g4
stw %xg4,[%sp+2259]
stw %xg5,[%sp+2263]
stw %xg6,[%sp+2267]
stw %xg7,[%sp+2247]
stw %xg8,[%sp+2251]
stw %g1,[%sp+2239]
stw %g3,[%sp+2243]
stw %g2,[%o0]
stw %g4,[%o1]
sub %sp,-240,%sp
retlなお、わかりやすいように遅延スロットの場所を変えるなどしている。
このアセンブリを読み解くと、関数のインライン展開は行われたが、関数ローカルに(すなわちスタックに)
作られた構造体が、後に不要だとわかっても削除できなかったことが推定される。
ワークアラウンドとしては、add関数の引数を参照渡しに変えたり、構造体にコンストラクタを書くことで
最適化が最後まで通り、富士通コンパイラでも即値で返せるようになる。
しかし、それを知らないユーザから見ると、自分の開発環境では全く問題にならなかった場所に余計なコードが大量に生成されており、知らないうちに性能が劣化している、ということになる。この現象をコンパイルメッセージから読み解くのは難しく、性能劣化に気づいた後にアセンブリを詳細に読み解かなければ、何が起きているかを把握できないであろう。
これについては、引数で受けた構造体に関する最適化処理が終盤で行われているため、そのステージでは ソースコードの情報が欠落しており、最適化が十分に適用されなかった旨が報告された。
これはC++言語に限らないことだが、大きなプログラムは一般に分割コンパイルが行われる。 分割コンパイルとは、更新のあったソースコードのみオブジェクトファイルを更新し、既にコンパイル済みの他のオブジェクトファイルとリンクをして実行ファイルを作成することで、全体のビルド時間の短縮をはかるものである。 数値計算、特にHPC分野で使われるコードは大規模化が進んでおり、フルビルドに何時間もかかるものが珍しくない。 したがって、分割コンパイルを行うわけだが、この時に問題となるのがプロシージャ間の最適化である。簡単な例を挙げよう。
整数を受け取って、その絶対値を返す関数adjustをadjust.cppに実装する。
int adjust(int a) {
if (a < 0)return -a;
else return a;
}他のファイルから使えるように、関数宣言をadjust.hppに書いておく。
#pragma once
int adjust(int a);これをmain関数から以下のように呼んで見よう。ファイルはmain.cppとする。
#include <cstdio>
#include "adjust.hpp"
int main() {
printf("%d\n", adjust(-5));
}さて、このコードをコンパイルする時にはadjustがどんな関数かわからないので、ただcallするしかない。g++でコンパイルした際、対応するアセンブリはこうなる(関数名が見やすいようにc++filtをかけてある)。
g++ -O3 main.cpp adjust.cpp movl $-5, %edi
call adjust(int)
leaq lC0(%rip), %rdi
movl %eax, %esi
xorl %eax, %eax
call _printfしかし、このコードを-fltoオプションをつけてコンパイルするとg++は即値を返せるようになる。
$ g++ -O3 -flto main.cpp adjust.cpp
$ objdump -S ./a.out
(snip)
0000000000400500 <main>:
400500: 48 83 ec 08 sub $0x8,%rsp
400504: be 05 00 00 00 mov $0x5,%esi
400509: bf 74 07 40 00 mov $0x400774,%edi
40050e: 31 c0 xor %eax,%eax
400510: e8 cb ff ff ff callq 4004e0 <printf@plt>
400515: 31 c0 xor %eax,%eax
400517: 48 83 c4 08 add $0x8,%rsp
40051b: c3 retq
(snip)callによる関数呼び出しが消え、adjust(-5)の実行結果である5がそのままprintfに渡されているのがわかる。
これは、リンク時最適化(Link time optimization, LTO)と呼ばれ、オブジェクトファイルを生成する際に中間情報を保存しておき、リンク時にあらためて最適化をかけてから実行バイナリを作るものである。
上記のコードは、たとえば分子動力学法において周期境界条件を補正するためのコードを模したものである。 最近のコンパイラは賢いこともあり、C++プログラマは、「この程度ならインライン展開されるだろう」と期待して、小さな関数をたくさん作ることがある。しかし、コンパイラがLTOに対応していないと、これらは全て関数コールになってしまう。特にループの中で何回も呼び出される関数がインライン展開されないと大きな実行ペナルティとなる。
LTOに対応していないコンパイラを用いる場合の簡単なワークアラウンドとしては、関数の実体をヘッダに書いてしまう、という方法が挙げられる。先程の例であればadjust.hppの中に関数定義を書いてしまえば富士通コンパイラも即値を返せるようになる。しかし、なんでもかんでもヘッダファイルに入れてしまうと分割コンパイルの意味がなくなってしまう。
また、いちいち「この関数はたくさん呼ばれるからヘッダファイルに入れよう」などと考えるのも面倒である。
本WGでの発表時点では富士通コンパイラはLTOに事実上未対応であった。実際にはCコンパイラはLTOオプションがあったのだが、あまり有効に機能していなかった。渡辺からLTOに対応して欲しいと要望を伝えたところ、千葉よりLTOが有効に機能するベンチマークがあること、コードの中小構文木の情報をオブジェクトファイルに埋め込むことで、ポスト京ではLTOに対応する予定である旨の報告があった。
渡辺が開発した分子動力学法コードは、並列化に疑似Flat-MPI法という手法を採用している。これはFlat-MPIの計算とMPI/OpenMPハイブリッド計算が、通信部分を除いて全く同じ内容になるという特徴を持っている。したがって、通信時間が無視できる条件においては、Flat-MPIとハイブリッド計算は同じ計算資源を使うならば同じ結果を与え、かつ計算時間もほとんど変わらないはずであった。事実、別のスパコンサイトでは両者に差は見られなかったが、「京」においてはシングルノードで実行した場合でもすべてをプロセス並列した場合と、スレッド並列した場合で、有意にマルチスレッド実行の方が遅い現象に悩まされていた。その後の調査により、スレッド並列で実行される関数内部で使われているstd::vectorが性能劣化の要因であることがわかった。問題となる関数を簡略化して抜き出すと以下のような処理となる。
#pragma omp parallel for
for(int i=0;i < thread_num; i++){
func(i);
}上記のように、スレッドごとに異なる引数で関数が呼ばれている。この関数の冒頭で、以下のようにstd::vectorが宣言され、利用されている。
void func(int thread_index){
std::vector<int> v;
//...
}分子動力学法のホットスポットは力の計算であるが、このケースも90%以上が力の計算に費やされており、上記の場所は別のサイトであれば全体の計算時間に比べて無視できるほどのコストしかかからないはずである。しかし、「京」でマルチスレッド実行すると上記の処理が極端に遅くなり、全体で10%程度の性能劣化が見られた。
しかし、この箇所を、
void func(int thread_index){
thread_local std::vector<int> v;
v.clear();
//...
}と、スレッドローカル宣言することで性能が改善し、Flat-MPIの場合とほぼ同じ実行時間となった。
C++言語を利用するプログラマのほとんどは、STL (Standard Template Library)を用いる。中でも、C言語の配列を置き換えるstd::vectorの利用頻度は非常に高い。std::vectorは、これまでプログラマが陽に行ってきたメモリ管理を隠蔽し、メモリの確保や解放を自動で行うことができる。その際、内部ではmallocやfreeといったメモリ確保、解放の関数を呼び出している。std::vectorが関数のローカル変数として使われていると、毎回malloc/freeが呼び出されることになる。しかし、thread_local宣言をすることでstd::vectorは暗黙にstatic修飾され、今回の実行条件ではfreeはプログラムの最後に一度だけ呼ばれるのみとなる。マルチスレッド実行では共有メモリとなることから、Flat-MPI実行に比べ、mallocはより大きなメモリプールについて管理を行う必要がある。これがマルチスレッド実行時の性能劣化原因であろうと推定される。
上記の問題の他、高並列時、およそ数百ノード規模の計算を行った場合、Flat-MPI実行に比べてハイブリッド並列実行が20%程度性能劣化すること、関連して、富士通独自の環境変数の役割が分かりづらいことについて報告した。性能劣化については完全な原因究明にはいたっていないが、環境変数の設定変更による性能変化の原因についてや、富士通独自のmallocの環境変数については富士通小田和が報告書をまとめているので参考にされたい。
本稿では、主に富士通C++コンパイラへの要望と、それに対する富士通の対応予定について述べた。本WGで指摘された問題点のほとんどは、ポスト京において改善される見通しである。よく言われることであるが、最適化について「銀の弾丸」は存在しない。一般論としては「環境を変えてもコンパイラがちゃんと最適化してくれるのでコードの可搬性が保たれる」というのは幻想に過ぎず、コードのホットスポットはハードウェアに応じてほぼ書き直す必要があるのが現状である。しかし、だからといってコンパイラの最適化機能はおざなりにしてはならない。むしろ、プログラマがホットスポットのチューニングに注力するためには、コンパイラが「露払い」として、それ以外の部分をきっちり最適化できなければならない。ポスト京を取り巻く環境を概観すると、コアの増加、NUMA、SIMD幅の増加といったハードウェアの複雑化に伴ってプログラマの負担は増え、プログラミング言語C++もC++11、14、17と機能が追加されていき、コンパイラ開発陣の負担も大きくなっている。プログラマも「コンパイラの最適化が甘い」と愚痴るのみならず、出力されたアセンブリを確認して「どんなコードを吐いて欲しいか」を開発に伝え、開発もそのフィードバックを受けてコンパイラを改善することで、ユーザと開発が両輪となり、より良いHPC環境を構築していくことが望ましい。
[1] "SIMD vectorization for the Lennard-Jones potential with AVX2 and AVX-512 instructions", H. Watanabe and K. M. Nakagawa, Comput. Phys. Commun., Vol. 237, pp. 1 (2019)
[2] https://github.com/kaityo256/mdacp
[3] メニーコア時代のアプリ性能検討WG 成果報告書 公開URL
http://www.ssken.gr.jp/MAINSITE/download/wg_report/mcap/index.html
本報告書に関連する発表スライドをSpeakerDeckにアップロードしてあります。