将来どうなるかはわからないのだが、とりあえず2018年現在、まだスパコンを使って数値計算をしようとしたらFortranかC++を使わなければならない。で、一般論としてコンパイラはFortranの方が最適化しやすく、似たようなコードをFortranとC++で書いたら、Fortranの方が早くなる、ということはまま起きる。それじゃ全部Fortranでかけばよいのか？というと、世の中そんなに単純では無い気がするよ、という感じのポエムです。
なお、以下に出てくる数値計算結果は筆者の(だいぶ昔の)経験を脚色したもので、実際に存在するアーキテクチャやコンパイラとは全く関係ありませんのであしからず。
分子動力学法(Molecular Dynamics method, MD)、特に短距離古典MDは、数値計算対象としてみると、二つの大きな特徴がある。一つは論理的に隣接するデータが固定していないこと。もう一つはループの内部で間接参照を伴うこと。もちろんこの二つは関連している。MDでは粒子が動き回るため、さっき相互作用してた粒子が遠くにいって、今は別の粒子と相互作用している、ということが起きる。空間を構造格子に切って計算する方法では、相互作用相手がずっと固定されているのとは大きな違いとなる(メッシュの切り直しとかはまた別の話ね)。また、相互作用相手が一定しないため、相互作用リストを作る必要がある。すると、相互作用を計算するためには、まずリストから相手の粒子番号を引っ張ってきて、その粒子番号で座標や運動量を取ってくる、という手続きを踏む必要があり、本質的にループが間接参照を伴う。これらはコンパイラによる最適化を強く阻害する。
で、MDは本当はそういう性質があるんだけれど、その二つをとっぱらったサンプルコードを書いてみる。つまり、全粒子対について相互作用を計算する。すると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とFortranは多次元配列の順序が異なるのでそれは入れ替えてある。さて、これをそのまま同じ石で、同じようなコンパイルオプションでコンパイル、実行すると、先程はかからなかった最適化がかかり、20%程度速度が向上する。
さて、Cで書いたコードを「そのまま」Fortranにするだけで、コードが数十パーセント早くなった。特にC/C++の特別な機能を使っていないにもかかわらず、である。「それじゃ最初から全部Fortranで書けばいいんじゃ？」と言いたくなるのをグッとこらえてもらって、ちょっと最適化してみよう。力の計算で、内側のループで値が変わらない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;
}
}すると、先程C++版にはかからなくてFortran版にはかかっていた最適化がかかり、さらに別の総和演算に関する最適化もかかって、実行時間は元のコードに比べて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同じ最適化をしているのだからあたりまえだが、同じ最適化がかかって同じ実行速度になる。Fortranの場合はもともと少し早かったので、この最適化による速度向上率は26%と、C版よりは向上率が低い。
実際に「全く同じコードをCからFortranに書き直すだけでかなり高速化される」という事実はある。しかし、多くの場合それではチューニングは不完全であり、ちゃんとチューニングすればどちらも同じ速度になる(同じようなアセンブリを吐くようにチューニングするから当然なのだが)。どこまでチューニングするかはもちろん場合によるが、そもそも「チューニングしよう」と思うからには「この最適化がかかるはずなのにかかってない/このくらいの性能が出るはずなのに出ていない」という知識があるはずで、その知識があればその最適化がかかるように/その性能がでるように最適化をするのはさほど難しいことではないはずで、「組み込み関数バリバリでアセンブリガリガリ」とか言い出さなければ、C++だろうがFortranだろうが、ちゃんと最適化すれば同じ石で同じような性能が出るはずなんですよ。
なので、「C/C++コンパイラがアホならFortranで書けばいいじゃん」というのは一種の思考停止だと思う。きちんとチューニングされたコードなら、(理想的には)言語によらず一定の性能が出るはずなんだから(まぁそうならない時もあるんだけどさ)。
ただし、「基本的にロジックだけ書いて、最低限の最適化しかしない」という立場を取ることも考えられる。その場合は(概ね同じコードを書いた場合にC++よりも最適化がかかりやすい)Fortranを使う、というポリシーは一定の合理性がある。
勘違いしないで欲しいんだけど、僕はFortranに比べてC++が優れているとは思っていない。例えばべき乗や複素数型が言語仕様にデフォルトで含まれていないなど、C++は数値計算言語として不満な点もある。僕がいいたいのは「C/C++コンパイラがアホならFortranで書け」というのはちょっと違うんじゃないかな、ということ。
ちなみに、今回のケースでは、最初のコードはC++用の特別なオプションを用いることでちゃんと最適化がかかって、Fortranと同等な性能が出る。でもまぁ、その後最適化するとさらに26%速度が向上するわけで。この26%という数字をどう思うかだけど、普段遣いのPCクラスタならどうでもいいとは思う。でもある程度大きな計算資源で並列化してぶん回す、というときに、数十パーセントの性能向上の余地のあるコードを投げる、というのはどうかな……
現状、スパコンの開発言語としてFortranかC/C++の二択、というのはそうなんだけど、最近はPythonやらJuliaやらを使ってもう少し幸せにコードがかける環境が整いつつある気がする。っていうか「下から上まで同じ言語で開発する」というスタイルは、もうプログラマのしんどさの限界を超えてるんじゃないかと思う。だいたい数値計算の重いところってある程度固定されているので、それをNumPy、もっと後ろはBLASみたいなライブラリ側で頑張って、上流は「組み上げ」に徹する、というのが幸せな気がするし、そういう未来は近づいて来ているんじゃないかな。ただ、最下層のチューニングはそれで良いとして、並列化はどうするんだろう？よくわかりません。
ちなみに、チューニングの話をすると必ず「そのうちコンパイラが賢くなるのでチューニングはオワコン」みたいなことを言い出す人がいる。僕の記憶にある限り、この言葉を15年以上前から聞いているのだが、それから15年たった今でも、そんな素敵な未来はあんまり見えてきてない。まぁでも、僕だってそういう未来が来てほしいし、そのために、つまり「いまあるコンパイラ」を、そういう「素敵なコンパイラ」に近づけるために、今日もコンパイラが吐いたアセンブリを見て「これはおかしいんじゃないの？」とかベンダーに報告を上げてるわけです。
あ、しつこいですけど、この文章はフィクションであり、現実のいかなるアーキテクチャやコンパイラ、人物とは無関係ですのでよろしくおねがいします。