3.4.4 科学計算におけるC++の状況
田中 義一委員
1．はじめに
C++には、科学計算に魅力的な機能がいくつもある。テンプレートを使った汎化プログラミング、記述性の高い演算子の多重定義、コード再利用のためのオブジェクト指向的機能である。しかし、この魅力的機能にもかかわらず性能に大きな問題があるため、科学計算の一般ユーザは、C++を本格的には使用していないように思われる[1]。実際に、C++を使用してプログラムを書くと、FORTRANに比べ何十倍も遅くなることはよくあることである。この状況は、命令レベル並列方式によって性能をあげるハイエンドプロセッサでは特にひどくなってきている。
本稿では、大きな抽象化ギャップに対するコンパイラによる最適化の困難さと、これを打破するためのアクティブライブラリ構築の試みを、具体的な例に基づいて紹介する。
2．変数のエリアス関係の増加
一般に、C++で記述すると、生データやポインタをクラスでラップすることが多くなる。例えば、以下のプログラムのように、CやFORTRANではデータdataaなどのポインタを直接アクセスしていたが、C++では小さな管理オブジェクトaを経由してアクセスすることが多い。
class A{double *dataa;}; class B{double *datab;}; class C{double *datac; };
void foo( A &a, B &b, C &c, int n){
for(int i=0;i<n;i++)a.dataa(i)=b.datab(i)+c.datac(i);}
性能のキーとなるループ内のロード・ストア命令の数を予想すると、最適化されたコードでは浮動小数点ロード2と浮動小数点ストア1つである。ところが、実際のC++コンパイラは、dataa,
datab,datacフィールド間にエリアスの関係にあることから、書き換えの可能性があると判断する。従って、ループ内にはこれらのフィールドをロードする命令が追加され、ロード命令は5つとストア命令が1つとなる。データを直接アクセスするのであれば｛
void foo( double *dataa,..) ｝、dataa,datab,datacは別名の関係にないため、そのような問題は発生しない。このように、C++で書かれたプログラムは、変数間のエリアス関係を形成しやすいので、冗長なコードや、並列度の低いコードが生成される。
エリアスの問題に対して、最近のコンパイラはコンパイルオプションや言語拡張で対応している。コンパイルオプションの例としては、エリアスに関する条件を緩和するオプションがあげられる。例えば、ポインタの自己参照がないとか、型の異なるポインタ同士は別名の関係にないなどである。すなわち、トリッキーなプログラムでなければ通常は守られているエリアス条件を仮定するコンパイラオプションである[2]。
ここでは、命令数の観点から述べたが、スーパスカラプロセッサやVLIWプロセッサ向けに並列度の高いコードを生成するためには、さらにポインタの指す先がエリアスの関係がないことが必要である。このために、restrictポインタと呼ばれるポインタ仕様がデファクトスタンダードとなりつつある。
3．大量の一時オブジェクト
C++では、一時オブジェクトの導入がたびたび必要になる。特に、多重定義演算子を含む式を処理する場合に発生しやすい。一時オブジェクトが導入されると、必要なコンストラクタやデストラクタ呼び出しや、オブジェクト間のコピー演算（関数の引数や返り値渡す）が発生するので性能上問題が多い。インライン展開により、一部のコピーは消去できることもあるが、複雑にネストされたオブジェクト間のエリアスの問題で削除できないものも存在する。コピーコンストラクタが定義されている場合に、名前付返却値（NRV）最適化などにより、関数の返り値に対してコピーの発生が抑止できることがある[3]。
以下に、一時オブジェクトの問題の典型として、STL（Standard Template Library）[4]の開発者であるStepanovが作成したベンチマークプログラムの説明と性能結果を示す。プログラムの内容は、「2000個のdouble型データの総和という単純な処理」をSTL風な記述をしたものである。このベンチマークの中には、大きなデータを扱うプログラムにおける、小さなオブジェクトの典型的な2つの扱い方が書かれている。すなわち、たくさんの小さなオブジェクトから構成される大きなデータセット（例えばdouble型配列）と、その大きなデータセットのアクセスするための小さなオブジェクト（例えばコンテナとイテレータ）の振る舞いである。ベンチマークは、以下に示す13個のループからなる[6]。
test0:ｃスタイルのループ
test1-test12: 関数オブジェクトを使用したSTLスタイル記述
template<class Iterator, class Number>
Number accumulate(Iterator first, Iterator last, Number result){
while(first!=last)result = plus(result, *first++); return result;}
test1,3,5,---- データは直接見えるdouble型変数 パラメータNumber <-double
test2,4,6,---- データはクラスDoubleの中のdouble型変数（ラッピング）
パラメータNumber <- Double
test1: 反復子として通常ポインタパラメータIterator <- *double
test2: 反復子として通常ポインタパラメータIterator <- *Double
test3: 反復子としてクラスdouble_pointer型変数（中にポインタ*doubleを含む）
パラメータIterator <- double_pointer
test4: 反復子としてクラスDouble_pointer型変数（中にポインタ*Doubleを含む）
パラメータIterator <- Double_pointer
．．．．．とtest12までつづく（図1参照）。
図１ test4のデータ構造とtest4のインライン展開後の擬似プログラム
accumulate(Double_pointer first, Double_pointer last, Double
result)
| while(oprator!=(&first, &last)){
仮引数をp1,p2と書く
| | int t1=0
| | t1 = ! operator=(p1, p2)
仮引数をp21,p22と書く
| | | return (*p21.current==*p22.current)
currentはDouble*
| | return t1
| Double_pointer t3;
| 次のoperator()関数の第３引数を求める関数;
| operator*( t3=operator++(int)(&first) ,&t3) operator*の仮引数をp5と書く
| | | Doube_pointer t4; oprator++関数の仮引数をp3と書く
| | | t4 =*p3;
| | | operator++()(&*p3) 仮引数をp4と書く 前置型＋＋
| | | | ++*p4.current
Double_pointer変数firstのcurrent
| | | | return p4
フィールドの++
| | | return t4;
| | return (*p5).current
*p5 ++前のDouble_pointerのコピー
| result=operator()(&plus, &result, 第３引数 Double*) 第２、第３仮引数をp7,p8
| ｜return operator+(p7,p8) 仮引数をp9,p10
| ｜ | double t10=0;
| ｜ | Double t11;
| ｜ | return Double( &t11,(t10=*p9.value+*p10.value, &t10)),t11
| ｜ | | value=*p11
Double コンストラクタの引数 double *p11
| ｜ | |
p11は加算結果
| }
| return result;
表1に、ある高性能スーパスカラプロセッサでの、３つの異なるコンパイラによるtest0-test4の性能と、ループ1回あたりの実行命令数を示す。コンパイラAは、C++の最適化では世界的に著名なコンパイラである。そのコンパイラでもtest0以外はあまりよい性能とはなっていない。一時オブジェクトは完全に除去されているため、無駄なコードはない。しかし、最近の計算機では、高速性の前提である命令レベル並列方式を生かすソフトウェアパイプラインに基づくコードが生成できなければ、命令数以上の性能差が発生する（本例は、総和と言う依存関係にある処理であるため、その傾向がさらに強くなっている）。このことは、同一周波数の他のスーパスカラプロセッサでも確かめることができた。メーカ製コンパイラを用いた場合、test0は714mflopsでループ1回あたり2.1個の命令、test4は17mflopsで7.5個の命令であった、命令数が3.5倍に増加しただけで、性能は1/42倍というひどい状態になっている。このように抽象度が高いソースプログラムに対して、合格点の最適化コードを生成するには、1つの見逃しも許されないということを意味している。
表1 Stepanovベンチマークの性能
|
|
コンパイラＡ |
コンパイラＢ |
コンパイラＣ |
|||
|
性能 |
命令数 |
性能 |
命令数 |
性能 |
命令数 |
|
|
test0 |
556mflops |
1.8 |
625mflops |
1.8 |
31mflops |
6 |
|
test1 |
91mflops |
5 |
556mflops |
1.8 |
33mflops |
6 |
|
test2 |
91mflops |
5 |
455mflops |
2.3 |
33mflops |
6 |
|
test3 |
89mflops |
5 |
21mflops |
10 |
33mflops |
7 |
|
test4 |
93mflops |
5 |
11mflops |
15 |
33mflops |
7 |
4．オブジェクト指向モデルと数値計算
C++で書かれた物理シミュレーションプログラムを見ると、以下のようなスタイルで書かれたプログラムを散見することが多い。このプログラム断片は、流体力学のシミュレーションにおいて単純なSOR法によって圧力を求める部分を表している。ここでは、2次元セル（i,
j）点での残差を求めており、この外側に領域を巡るループがある。
Poisson::zansa( int i, int j)
{ double diffx = 0, diffy=0;
if( grid().cell(i,j).isInsideFluid() )
{ diffx = grid().P(i-1,j ) - 2 * grid().P(i,j) + grid().P(i+1,j
);
diffy = grid().P(i ,j-1) - 2 * grid().P(i,j) + grid().P(i
,j+1); }
else {
double p = grid().P(i,j);
if( grid().cell(i,j).isFluid(Cell::DOWN )) diffy += grid().P(i
,j-1) - p;
if( grid().cell(i,j).isFluid(Cell::LEFT )) diffx += grid().P(i-1,j
) - p;
．．．．．}
return diffx * dx + diffy * dy - so()(i,j);
}
この断片から、セルi,jの状態は、すべてi,jの物理量を表すオブジェクトが知っているというモデルで書かれていることがわかる。従って、各セルごとに、「このセルは流体を含んでいますかisInsideFluid」、「このセルの左隣は流体を含んでいますかisFluid(Cell::LEFT)」、「隣のセルの圧力はＰ(i+1,j)」、、、というメンバ関数による問い合わせが発生することになる。このようなプログラムの機械語レベルの構造を見ると、最内側ループは、物理量の属性（計算のための属性）の種類だけ分岐処理を含んだ複雑な処理になる。
スーパコンピュータ向けにプログラムをFORTRANで書く場合、高速化を望むプログラマは、「最内側ループはできるだけ単純」に書くというプログラム指針に従っていた。単純という意味は、次のループ繰り返しi+1での処理が、データを読まなくても静的に予測（線形）[3]できるということにつきる。上記の例題が長方形の領域の場合には、内部領域内と特殊処理の境界を別々のループにして、内部領域内は性能向上が達成できるシンプルなループ構造にするアプローチがとられていた。すなわち、内側ループにある条件文を外側に追い出したことに相当する。また、任意の形状の場合であっても、間接ベクトルを使用して、同様な処理を記述することが行われてきた。逆に、FORTRANでは、徹底的にセルi,
jに着目したオブジェクト指向モデルで書くのは煩雑すぎると思われる。これが、FORTRANがスーパコンピュータと相性がよい理由と考えられる。
上記の例題のような物理問題において、解法に対して素直なオブジェクト指向モデルによって設計されたプログラムを、FORTRANのごとく書くのは、データと手続きが融合されているので困難と思われる。高速化のためには、物理解法モデルだけでなく、計算の規則性に基づいて計算領域を分割するような、計算の特性に基づくクラス設計を導入する必要がある。
5．式テンプレート
ベクトル、行列、配列演算が多くの科学計算の基礎である。C++では演算子多重定義により、各クラスの演算をカスタマイズできる。これにより、STL[4]などでは
void foo(vector<complex> &y, vector<complex> &x, complex
a) {
y += a*x; }
のように、積和ベクトル演算を非常に自然な方法で表現できる。しかし、演算子の多重定義の評価は2項演算であるので、上記の式では、まずt1= a*x、次にt2=y+t1を行い、最後にy=t2と代入される。ここで、t1,t2のオブジェクトは中間結果のための一時オブジェクト（ベクトル）である。2項演算の式評価では、1つのループで式を評価するのではなく、3つの分かれたループで評価するためにオーバヘッドが大きい。また、一時オブジェクトは動的に取られるために、さらに性能が悪化する。例えば、上記スタイルと、Cによる低レベルの記述による性能は、ベクトル長が1000の場合、あるスーパスカラプロセッサでは、「ベクトルクラスで記述333Mflops」「低レベルで記述640Mflops」であった。
上記の問題は、FORTRANのアレイ演算でも発生することがある。しかし、通常コンパイラが、FORTRANからローレベルの内部言語に変換する際に、左辺と右辺でエリアスの可能性がないと判断した時は、1つのループにすることが一般的である。同様のことをC++で行うことは一般的に困難と考えられる。FORTRANのアレイ演算は、言語の組み込みであるのに対して、C++ではユーザ（クラスライブラリ）が定義するものであるからである。コンパイラのフロントエンドが、各演算子がどのメンバ関数での定義であるかを見つけて、必要なら一時オブジェクトを生成しながら対応する関数を生成する。従って、C++の場合は、一度複数ループを生成してから、再度インライン展開後に隣接するループを融合する方法を取らざるを得ない。しかし、ループの間には、例外を投げうる一時オブジェクトの動的割り当て文の存在もあり、ループ融合ができるかどうかを判定する一般的方式は考えにくい。
上記の問題に対して、blitzクラスライブラリ[5]では、2項ごとの評価を式テンプレートと呼ばれる技術で解決した。以下では、簡単にdouble型のベクトル加算演算だけが可能なクラスライブラリのみを示す。ここでは、明らかなコンストラクタやデストラクタ等や、public/privateは削除してある[4]。
template<class Op>
class vector{
vector & operator = (Vop<Op> const &expr){
(1)
for(int k = 0; k != s_; ++k ) a_[k] = expr[k]; return *this;
}
double const& operator[](int k) const { return a_[k]; }
double& operator[](int k) { return a_[k];}
double *a_;
int s_;
};
template<class Op>
class Vop {
double operator[](int k) const { return op_[k]; }
int size() const { return op_.size();}
Op ope() const { return op_; }
Op const op_;
};
template<class Va, class Vb>
class Vadd{
double operator[](int k) const {return a_[k] + b_[k];}
int size() const {return a_.size();}
Va const &a_; Vb const &b_;
};
Vop<Vadd<vector,vector > > operator+(vector const &a,
vector const &b){
(2)
return Vop<Vadd<vector, vector> > (Vadd<vector,vector
>(a,b));
}
template<class Va>
Vop<Vadd<Va,vector > > operator+(Vop<Va> const &a,
vector const &b){
(3)
return Vop<Vadd<Va, vector> > (Vadd<Va,vector>(a.ope(),b));
}
上記のクラスライブラリが与えられたとき
vector x(1000),y(1000),z(1000),w(1000)
x=y+z+w:
において、ベクトル演算が効率的に実装されることを示す。まず、y+zを考える。y,zがvector型なので、(2)のoperator+が使われ、Vadd<vector,vector>オブジェクトを格納するVop<Vadd<vector,vector> >型のオブジェクトが作られる。ここで、この加算はVop型の小さなオブジェクトを作ることだけを行う。y+zの評価は先送り（遅延評価）にされている。このオブジェクトは閉包（closure）と呼ばれるものと密接に関係している。同様に (y+z)+wはVop<Vadd<vector,vector> >型とvector型の加算なので(3)のoperator+が使われ、Vop<Vadd<Vadd<vector,vector>,vector> >型オブジェクトを生成する。すなわち、以下のようにコンパイル時にテンプレートの具現化が行われる。
=y+z+w;
= Vop<Vadd<vector,vector> >( Vadd<vector,vector>(y,z))+w
（yz= Vop<Vadd<vector,vector> >( Vadd<vector,vector>(y,z))とすると）
= Vop< Vadd<Vadd<vector,vector>,vector> >( Vadd<Vadd<vector,vector>,vector>
(yz.op_ , w)) （=yzwとおく）
そして、これがvectorに代入されたときに、(1)の代入演算により添え字演算子が以下のように展開され、演算評価が行われる（ここまで評価が遅延された）。このため、従来のvectorクラスを使用した時のように、ループが分割されることはなく、Cにより低レベルで記述した場合と同様な性能が実現できることになる。
x.operator=(Vop<Vadd<Vadd<vector,vector>,vector> > const
&expr
) for(int k = 0; k != s_; ++k) a_[k] = expr[k]
= yzw.op_[k] = yz.op_[k]+w[k]
=y[k]+z[k]+w[k]
6．アクティブライブラリ
ここまで述べてきたように、C++で書かれた抽象化のレイヤを、従来のようなコンパイラの最適化（transformational approach）のみで解決することは難しい。コンパイラはプログラムの意味がわからないので、できることが限られるためである。5．で述べたblitzライブラリ[5]は、抽象化レイヤに対する最適化技術をライブラリで行うことを目的としている。ここでライブラリとは、従来の意味での種種の機能の集積である伝統的なパッシブなライブラリでなく、C++の持つ高度のテンプレート機能により、コンパイル時に、コードをオプティマイズするアクティブなライブラリである。すなわち、最適化コードはテンプレート技術を使用してライブラリ自身によって生成（generative
approach）される。
5.で述べた式テンプレート技術は、評価をする前にいくつかの情報を1つの関数にまとめなければならない性質をもつ様様な問題にて応用ができる。
例えば、
Array<double, 3> A,B,C,D;
A=B+C+D;
とすると、まず、5.で述べたように、式テンプレートが式のパースツリーを生成することで、テンポラリ変数を生成させないで以下のように3重ネストループを生成できる。
for(int i=0; i< N1; i++)
for(int j=0;j<N2;j++)
for(int k=0;k<N3;k++)
A(i,j,k)=B(i,j,k)+C(i,j,k)+D(i,j,k);
そしてパースツリー操作の際に、blitzライブラリでは、次のような最適化を行っている。ループ交換、ループ一重化、アンローリング、配列ステンシル（5点差分、7点差分など）を検出したあとのタイリング[5]などである。アレイクラスにおいてインデクスをトラバースするイテレータをいくつか用意して、適用しているように思われる。これらの最適化は、高レベルに抽象化されたソースプログラムに適用されるので、汎用性が高い。一方、これらの最適化は、コンパイラでも行われるが、低レベルに変換されたソースプログラムへの適用であるため、汎用性は小さい。
7．おわりに
C++をFORTRANなみの性能で使うにはまだ障害が多い。今後は、C++のような抽象度の高い言語に対する最適化は、コンパイラではなく、テンプレート機能によるアクティブライブラリによる最適化が有望である。ユーザは、これらの最適化されたライブラリを使って、初めて高速性を享受できる可能性がある。
blitzのようなライブラリが機能拡張[6]，標準化され、商品化されることが望まれる。なお、テンプレートを駆使したアクティブライブラリのソースプログラムを理解することは非常に困難で、一般のユーザが自分で構築することは無謀と思われる。また、テンプレートプログラミングの最大の問題は、膨大なテンプレートの具体化のため（部分計算）コンパイル時間が膨大になることである。何らかの方策（プリコンパイルヘッダなど）が考えられないと、実用的にはならないと思われる。
参考文献
[1] http://www.extreme.indiana.edu/hpc++/index.html
[2] http://ww.acl.langl.gov/pooma/
[3] S.B.Lippman：C++オブジェクトモデル，トッパン(1997)
[4] D.R.Musser, and A.Sani: STL Tutorial and Reference Guide, Addison-Wesley (1996).
[5] http://www.oonemerics.org/blitz
[6] http://www.kai.com/benchmarks/stepanov
[7] D. Vandevoorde: C++ Solutions Companion to the C++ Programming Language, Addison-Wesley (1998).
[1] 研究者レベルでは、C++の持つオブジェクト指向プログラミング機能を生かした並列アプリケーション記述の研究が盛んである。例えば、C++の言語仕様を変更することなくC++の持つテンプレートおよびクラス機能を用いて並列処理を実現するアプローチをとるHPC++、同じくテンプレートやクラス機能を用いてアプリケーション寄りのライブラリを提供するPOOMAがあげられる[1][2]。
[2]また、過去のトリッキーなプログラムの実行は保証せずに、ポインタに関して緩和したエリアス条件を前提としてコードを生成するコンパイラもある。
[3]長方形領域の「右隣のセル圧力の大きさはＰ(i+1,j)」ですら、ループ内では予測できない。領域は、計算上、一般に内部領域と境界領域に分けられるが、セルが境界点にある場合は境界条件モデルによって、右隣の定義が異なる。
[4]ここで示したコードは実際のblitzの実装とは異なる。式テンプレートをわかりやすく説明するためのものである[7]。
[5]空間のトラバースの例としてHilbert curveによるものがあり、これを配列ステンシルに適用すると、大規模プログラムにおいて、キャッシュの再利用性をすこし高める実験結果が得られている。もちろん、タイリングが望ましいが、2次元までしかまだ対応していないようである。
[6] blitzライブラリは疎行列や、対称行列には対応していないよう。