HPCプログラマの書くコードってなんで古臭いの？

TL;DR

Q: HPCプログラマの書くコードってなんで古臭いの？
A: そうしないと性能が出なかったから

はじめに

スパコンを使ってそれなりに長いこと研究をしています。初めてスパコンを使ったのは今から25年くらい前です。発展の早い世界ですので、昔と今ではだいぶ違ってきています。この前、スパコンを使う時の「常識」がすでに現在では常識ではなくなっているっぽいことを知り、私も「そっち側」になったのか、と感慨深くなりました。以下は、HPC業界の語り部に片足を突っ込んだおっさんの戯言です。

なお、以下では主に「京」の話をするため、富士通のC++コンパイラに文句を言う形になりますが、富士通が悪いのではなく、昔のスパコンに搭載されていたC++コンパイラはどれもかなり酷かったということはあらかじめ言っておきたい気がします。IBMのコンパイラとかも酷かったし、なんならIntelのコンパイラもバグだらけだったりしたし。

C++とFortran

Fortranというと、古臭いプログラミング言語の代表例としてよく紹介されます。「スパコンでは、まだFortranが現役で使われている」と聞くと、驚く人も多いようです。しかし、昔、全く同じコードを書いても、C++では性能が出ず、Fortranではそれなりの性能が出る、ということがありました。昔と言っても、「京」コンピュータの頃(2012～2019年)でもそういう状況だったので、個人的には「すごく昔」というわけではないのですが。

私の専門が分子動力学法なので、単純な分子動力学法コードを書いてみます。本来ならカットオフを導入して隣接粒子リストを作る必要があるのですが、重力N体のように$O(N^2)$で書いてみます。こんな感じです。

#include <cstdio>
const int N = 20000;
const int D = 3;
const int X = 0;
const int Y = 1;
const int Z = 2;
double q[N][D],p[N][D];
void
calcforce(void){
  for(int i=0;i<N-1;i++){
    for(int j=i+1;j<N;j++){
      double dx = q[j][X] - q[i][X];
      double dy = q[j][Y] - q[i][Y];
      double dz = q[j][Z] - q[i][Z];
      double df = (dx*dx + dy*dy + dz*dz);
      p[i][X] += df*dx;
      p[i][Y] += df*dy;
      p[i][Z] += df*dz;
      p[j][X] -= df*dx;
      p[j][Y] -= df*dy;
      p[j][Z] -= df*dz;
    }
  }
}
int
main(void){
  for(int i=0;i<N;i++){
    q[i][X] = i;
    q[i][Y] = i;
    q[i][Z] = i;
    p[i][X] = 0.0;
    p[i][Y] = 0.0;
    p[i][Z] = 0.0;
  }
  calcforce();
}

C++というか、ほとんど生のCなんだけれども、これをコンパイルするとソフトウェアパイプライニングが効かず、性能が出ません。「京」は豊富なレジスタを武器に、ループをソフトウェアパイプライニングで力任せに展開して演算のレイテンシを隠す、というアーキテクチャなので、コンパイラでソフトウェアパイプライニングができないと大きく性能が落ちてしまいます。

さて、同じコードをFortranで書いてみましょう。

subroutine calcforce(q,p)
implicit none
integer,parameter:: N=20000
integer,parameter:: D=3
integer,parameter:: X=1,Y=2,Z=3
double precision:: q(D,N), p(D,N)
integer :: i,j
double precision:: dx, dy, dz, df
do i=1,N-1
  do j=i+1,N
    dx = q(X,j) - q(X,i)
    dy = q(Y,j) - q(Y,i)
    dz = q(Z,j) - q(Z,i)
    df = dx*dx + dy*dy + dz*dz
    p(X,i) = p(X,i) + df*dx
    p(Y,i) = p(Y,i) + df*dy
    p(Z,i) = p(Z,i) + df*dz
    p(X,j) = p(X,j) - df*dx
    p(Y,j) = p(Y,j) - df*dy
    p(Z,j) = p(Z,j) - df*dz
  enddo
enddo
end
 
program force
integer,parameter:: N=20000
integer,parameter:: D=3
integer,parameter:: X=1,Y=2,Z=3
double precision:: q(D,N), p(D,N)
integer :: i
do i = 1, N
  p(X,i) = 0.0
  p(Y,i) = 0.0
  p(Z,i) = 0.0
  q(X,i) = i
  q(Y,i) = i
  q(Z,i) = i
enddo
call calcforce(q,p)
end program force

特に特別なことはしていません。計算ロジックは全く同じです。ただし、C/C++とFortranは多次元配列の順序が異なるので、そこだけ入れ替えてあります。これを「京」でデフォルトの最適化オプションでコンパイル、実行すると、先ほどのC++では効かなかったソフトウェアパイプライニングがかかり、速度が20%以上向上します。

スパコンにおける20%の性能向上というのはかなり大きいです。京では1000ノード24時間の計算とか普通に実行するため、20%性能向上すると、4800ノード時間だけ浮くことになります。何本もジョブを投げるので、それだけ性能向上はコスト削減に直結します。

さて、Cで書いたコードを、少しだけ手で最適化してやりましょう。力の計算で、内側のループで値が変わらないi粒子の座標を外に出し、かつi粒子の運動量を毎回書き戻さないように外に出してみます。

void
calcforce(void){
  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;
  }
}

すると、前回はかからなかったソフトウェアパイプライニングが適用され、性能がもとのコードに比べて40%近く向上します。この時点でFortranよりも早くなります。

しかし、Fortranでも同様な最適化をかけると、当たり前ですが同じくらいの速度になります。

subroutine calcforce(q,p)
implicit none
integer,parameter:: N=20000
integer,parameter:: D=3
integer,parameter:: X=1,Y=2,Z=3
double precision:: q(D,N), p(D,N)
integer :: i,j
double precision:: dx, dy, dz, df
double precision:: qx,qy,qz,px,py,pz
do i=1,N-1
  qx = q(X,i)
  qy = q(Y,i)
  qz = q(Z,i)
  px = p(X,i)
  py = p(Y,i)
  pz = p(Z,i)
  do j=i+1,N
    dx = q(X,j) - qx
    dy = q(Y,j) - qy
    dz = q(Z,j) - qz
    df = dx*dx + dy*dy + dz*dz
    px = px + df*dx
    py = py + df*dy
    pz = pz + df*dz
    p(X,j) = p(X,j) - df*dx
    p(Y,j) = p(Y,j) - df*dy
    p(Z,j) = p(Z,j) - df*dz
  enddo
  p(X,i) = px
  p(Y,i) = py
  p(Z,i) = pz
enddo
end

同じ最適化をしているのだから当たり前なのですが、コンパイラはほとんど同じアセンブリを吐き、結果として同じ実行速度になります。

さて、これを見てわかるように、スパコンでは全く同じ計算コードをC++で書くか、Fortranで書くかで大きく性能が変わることがよくあります。この例ではその差は20%程度でしたが、「京」以前のスパコンでは、そもそも全く勝負にならない、ということもよくありました。

もちろん、ここで示したように、最終的に同じアセンブリを吐くようにすれば同じ結果は出ます。しかし、先程のコードを見て、コンパイラメッセージから「ここが原因でソフトウェアパイプライニングできていない」と判断し、上記の様に修正するのは、よほど慣れたプログラマでないと難しいでしょう。そもそもホットスポットが一つだけならともかく、ホットスポットが分散しているようなコードで一つ一つ原因を探っていくのは面倒だし、それをやったからといって研究が進むわけではありません。「それなら最初からFortranで書いておけば楽じゃない？」と思うのも無理はないでしょう。

私は修士論文研究ではC++でコードを書きましたが、そのうちCになり、最終的にはFortranで書き直したコードで研究をするようになりました。途中からC/C++を使うようにしましたが、ずっとFortranを使ったほうが楽だ、という人の気持ちはよくわかります。

クラスの利用について

C++を使う時、機能やデータ構造をクラスにまとめる、ということをよくやると思います。しかし、それをすると変なところで最適化が阻害され、性能が出ないことがありました。

私のコードでは、分子動力学計算に必要な変数をVariablesクラスにまとめて使う、という設計になっていました。以下のように、座標と運動量、定数をまとめて、

class Variables{
  public:
    double C0;
    double q[N][D],p[N][D];
};

と定義しています。Nが原子数、Dが次元、座標がq、運動量がpです(ハミルトンダイナミクスを解いている気持ちなので、変数がx,vではなくq,pになっています)。

さて、ベンチマークコードでは性能が出るのに、本番コードでは思ったような性能が出ない、という問題に直面し、調査していったところ、コンパイラの最適化機能がクラス内のメンバ変数の定義順序に依存する ことがわかりました。先ほどの例と同様な、以下のようなコードを考えます。このコードでは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;
  }
}

これはVariablesクラスのインスタンスを受け取り、距離から定まる適当な力を計算して運動量に足し込むコードになっています。先ほど出したサンプルコードを、Variablesクラスのインスタンスを受け取るように変更しただけです。 冒頭で座標や運動量のポインタを定義して代入していますが、これは(少なくとも京稼働当初の富士通C++コンパイラでは)こうしないと最適化がかからなかっためです(闇)。 さて、このコードをそのままコンパイラに食わせるとソフトウェアパイプライニングがかからず、性能がでません。しかし、ここでVariableクラスの宣言順序を入れ替えてみます。

class Variables{
  public:
    double q[N][D],p[N][D];
    double C0; // 宣言を後ろにずらした
};

すると力の計算ループにソフトウェアパイプライニングが適用され、実行性能が倍以上向上しました。クラスの宣言順序を入れ替えただけで、です。今回のケースをユーザから見ると、「使われないクラスのメンバ変数の定義位置により、プログラムの性能が倍以上変わる」という奇妙な結果となっています。コンパイラ開発チームによると、こちらはrestrict修飾されたポインタ間の最適化処理に問題があったのが原因だったそうで、今は改善されているとのことです。

先ほどのコードもそうですが、普通にデバッグ or チューニングしていて「使われないクラスのメンバ変数の定義位置により、プログラムの性能が倍以上変わる」ことに気がつくのは難しいし、そもそも数値計算屋にそこまで突き詰めて調査させるべきか、という問題があります。最初からクラスを使わずに生配列で書いておけば全く問題はおきないわけで、このような事例を何度も目にしたスパコンプログラマが「なるべくクラスを使わない」ポリシーになったとしても、責められるべきではないでしょう。

STLの利用について

C++を利用する大きな動機の一つにSTLを使いたいから、というものがあります。生配列ではなくstd::vectorを使ってコードを書きたい、と思うのは普通でしょう。しかし、先程のコードを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にしても同様です。こちらはコンパイラチームによりSTLの内部で使われているvolatile変数が最適化を阻害していることが報告されました。こちらについても改善を検討中とのことで、おそらくLLVMベースの富岳では改善されているものと思います(すみません、未検証です)。

Ranged-base for

サイズがコンパイル時にわかっているなら、ranged-base forを使いたくなります。例えばstd::vector<int> vが保持する要素全てを2倍にするコードは

  for (unsigned int i = 0; i < v.size(); i++) {
  v[i] *= 2;
}

とも書けますが、C++11からなら

  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;
  }
}

ある配列の各要素について、絶対値を代入するだけのコードです。このコードをコンパイルすると、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命令が発行されていることが確認できます。

こういう関数一つを切り出せば、何が起きているかが予想できるので、対処も難しくないでしょう。しかし、これが数万行もあるコードの一部だったら？性能劣化の原因がranged-base forであることに気づけるでしょうか？

すごい苦労の上で、ranged-base forが原因であることがわかったプログラマは、次に何をするでしょうか？「なるべくC++11の機能を使わない」と学習してしまうのではないでしょうか？

マルチスレッド環境下でのstd::vector

私が京で最も苦労したのは、「全く同じコードがflat-MPIとハイブリッド実行で性能が大きく異なる」ことでした。分子動力学法においては、力の計算の最内ループをOpenMPでスレッド並列化し、さらにMPIを用いて領域分割することでハイブリッド並列とすることが多いです。しかし、私は「疑似flat-MPI法」という手法を開発し、スレッド並列においても領域分割を採用しました。これによって、flat-MPIの場合とハイブリッド実行の場合で、コアが担当する計算が全く同じになる、というメリットがあります。したがって、通信処理が無視できる領域では、flat-MPIとハイブリッドで全く同じ性能になるはずです。ところが、「京」ではシングルノード実行でも、8プロセスx1スレッド実行と、1プロセスx8スレッド実行に大きな差がつく、という問題がありました。

この問題のやっかいなところは、全体の実行時間は確かに遅くなっているのに、各関数ごとにストップウォッチを使ってウォールタイムを測ると、どこも遅くなっていないことです。プロファイルをとってもよくわからない関数が重くなっていたので、各関数に

void func(){
stopwatch s("func");
s.start();
//処理

s.stop();
}

のようなストップウォッチを追加し、時間を測定したのですが、これだとflat-MPI実行とハイブリッド実行で性能に差が出ません。

最終的に、これは「マルチスレッド実行時におけるstd::vectorのデストラクタが遅いのが原因」と判明しました。マルチスレッド実行時に各スレッドでstd::vectorがmallocを呼び、かつMPIの通信が混ざり、その後でスレッドがjoinしたタイミングでまたstd::vectorを使うと、シングルスレッド実行されたstd::vectorのデストラクタが極めて遅くなっていました。デストラクタが原因だったので、先程の関数内にストップウォッチを仕込んでもひっかかりません(デストラクタは関数終了時に呼ばれるから)。

わかるかっ！！！

当然、マルチスレッド環境下におけるstd::vectorの利用(より正確には、マルチスレッド環境下での頻繁なmalloc呼び出し)が問題だったので、生配列を使えばこういうことは起きません。

こういう事例を何度も経験したHPCプログラマが「なるべく生配列を使おう」というポリシーになったとして、責められるべきでしょうか？

まとめ

私が真面目に問題解決に取り組んだのが「京」だったため、富士通C++コンパイラへの愚痴のようになってしまいましたが、富士通C++コンパイラは他のベンダーのC++コンパイラに比べればかなり良い方で、昔はもっとひどい状態だったように思います。ただし、その時には私はすぐに諦めてFortranで書き直したり、生配列に書き直したりして対応していたため、何が原因かまでは調査していません。また、富岳ではLLVMベースになったことで、上記の問題のほとんどは解決しているものと思います(検証していないのでわかりませんが)。

私が言いたいことは、たった10年ほど前に、上記のようなことが起きていたということです。10年前は「最近」とは言えませんが、「大昔」とも言えないと思います。にもかかわらず、「なるべくクラス使わない」「なるべく生配列を使う」「なんならFortranで書き直す」といった対応策が「古臭い」と言われてしまうのは、ちょっと悲しいです。

現在のスパコンでは、このようなことはもう起きていないのかもしれません。しかし、20年近く「新しめの言語仕様を使うとひどい目にあう」ことを学習してきたHPCプログラマが、古臭いコードを書いたとして、それをバカにするのはやめてあげてください。HPCプログラマの心はそんなに強くありません。

あと、HPCプログラマは開発環境に文句を言うだけではダメで、ちゃんと「どういうアセンブリが出るべきなのか」とセットで開発部隊に報告書を上げるべきだと思います。ユーザとベンダは両輪となってHPC環境を良くしていかなければいけません。ユーザは「お客さん」ではないのです。

この小文を見て、私と同じような古参HPCプログラマが「うんうん」と頷いてくれることを信じて筆を置きます。

参考文献

• 富士通C++コンパイラの最適化機能の改善について
• 短距離古典分子動力学法のSIMD化及びC++コンパイラへの要望について
• 「京」におけるマルチスレッドmalloc

追記(2025年1月4日)

軽い気持ちで書いた記事ですが、思ったより反響がありました。タイトルの付け方があまり良くなかったかな、と反省しています。ここで勘違いして欲しくないのは、私はコンパイラには文句を言っていますが、コンパイラ開発部隊に文句を言いたいわけではないことです。上記の事象については全て富士通側に報告し、富士通側の回答も含めてサイエンティフィック・システム研究会の「メニーコア時代のアプリ性能検討ワーキンググループ」の成果報告書として公開されています。開発部隊がユーザの声にちゃんと真摯に答えている、ということは知っておいて欲しいと思います。

言語処理系を開発・運用するにはかなりのコストがかかります。一方、ハードウェアの「理論ピーク性能」であったり「電力性能」といった指標に比べて「言語処理系の成熟度」や「最適化性能」は具体的に数値化することが難しく、セールスポイントにすることが困難です。しかし、ユーザが「コンパイラはハードを買うと無料でついてくるもの」という意識でいると、ベンダーもコンパイラにコストをかけてくれません。処理系の開発では、言語の新仕様の追従やバグ修正が優先され、最適化については後回しにされがちですが、それでもユーザが声を挙げていかなければ状況は改善しません。コンパイラ開発部隊が上に人を増やせというには、溜まったチケットの数(もしくは対応を断念してクローズしたチケットの数)など、客観的な数字を見せるしかないからです。

現在、C++コンパイラをフロントエンドからバックエンドまでまるっと開発できるベンダーは数少ないと思われます。富岳はLLVMベースになり、それはそれで良い判断だと思いますが、全ての言語処理系がLLVMエコシステムに依存するのが健全な未来であるとは思えません。どこかにゼロからコンパイラ開発ができる「芽」を残して置かなければなりません。

記事本文にも書いたことですが、ユーザは「お客さん」ではありません。ベンダーと一緒に環境を良くしていくために一緒に努力しなければなりません。なにか処理系に不満があったら、単に遅いとか文句を言うのではなく、「どういう入力に対してどういう出力であるべきか」をセットにして報告を挙げなければなりません。多くの場合、開発部隊には余力がなく、対応が難しいかもしれません。しかし、声を挙げ続けなければ、状況は改善しません。

この記事が、単に「ベンダーの独自開発コンパイラはダメだよね」と捉えられるのではなく、「だから一緒に改善していかないとダメだよね」と認識して欲しいと思い、蛇足ながら追記しました。