uBLAS Overview

Rationale

あらゆる種類の数学的なソフトウェアが効率を失することなくC++ 言語で書かれ得るのなら、それは素晴らしいことだ。しかし C++ 型システムを犯すことなくこれを達成する何かが発見されない限り、 Fortran やアセンブラ、アーキテクチャ特定の拡張に頼る方がいいだろう (Bjarne Stroustrup)。

この C++ ライブラリは行列とベクトルに関する基本的な線形代数の構築、及びそれらに対応する抽象的な演算といった水準での科学技術計算を目指している。基本的な設計目標は、以下の通りである:

数学的表記
効率
機能性
互換性

もう一つの目的は、そのような行列とベクトルクラスを使うことによる抽象化の不利益が、許容できるものかどうかを評価することである。

Resources

このライブラリの開発はいくつかの似たような努力により導かれた:

BLAS by Jack Dongarra et al.
Blitz++ by Todd Veldhuizen
POOMA by Scott Haney et al.
MTL by Jeremy Siek et al.

BLAS は最も広く使われている、基本的な線形代数構築のためのライブラリだろう。そのため、これはデファクトスタンダードと呼ばれている。そのインタフェースは手続き志向であり、個々の関数は単純な線形代数演算からはある程度抽象化されている。 Fortran とその最適化を使って実装されてきたという事実により、これはまた、最も高速なライブラリのひとつだろう。我々はライブラリをオブジェクト指向で設計し実装することにしたので、技術的なアプローチは異なる。しかし、我々のライブラリの言葉で BLAS の抽象化を表現することや、これらの実装の効率を比較することは可能であるべきだろう。

Blitz++ は C++ で実装された印象的なライブラリである。その主な設計は、多次元配列と、テンソルを含むそれに関連する演算子を志向している。 Blitz++ の作者は、ライブラリは式テンプレートとテンプレートメタプログラミングを使った実装上のテクニックによって、同じような Fortran のコードと同等かそれ以上のパフォーマンスを達成した、と述べている。しかし、我々が、我々自身の設計と実装のアプローチで開発するのには多少理由がある。誰もが伝統的な線形代数や、他の数学的アルゴリズムを、 Blitz++ を使って実装しようとするかどうか、我々は知らない。また、現在でも Blitz++ は実装のイディオムにより、最も進んだ C++ コンパイラを必要とするのである。一方 Blitz++ は、抽象化の不利益を許容可能な程度におさめるために、式テンプレートの使用が不可欠である、ということを我々に解らせてくれた。

ROOMA の設計目標は多くの点で Blitz++ と似ているだろう。これは Blitz++ のクラスに関するコンセプトを、偏微分方程式と理論物理学の領域から拡張した。実装は同じようなアーキテクチャをサポートしている。

MTL は C++ で基本的な線形代数演算をサポートする為のもう一つのアプローチである。その設計は主に、 BLAS と C++ 標準テンプレートライブラリの影響を受けている。線形代数ライブラリは BLAS と同等の機能性を提供しなければならない、というその洞察については、我々も同じである。一方で我々は、 C++　標準テンプレートライブラリのコンセプトは、必要とされる数学計算をサポートするとは証明されていないと考えている。もう一つの違いとして、 MTL は現在、式テンプレートを使っていないように見える。これは 2 つの結果の内の 1 つをもたらす: つまり、表現力を失う可能性と、パフォーマンスを失う可能性である。

Concepts

Mathematical Notation

数学的表記の用法は科学技術アルゴリズムの開発を容易にするだろう。そこで、基本的な線形代数のコンセプトを実装している C++ ライブラリは、選択された C++ 演算子を行列やベクトルクラスに対して注意深くオーバーロードするべきである。

我々は以下の通り原初的なものに対して演算子オーバーロードを利用することにした。

説明	演算子
ベクトルと行列のインデックス	`vector::operator(size_t i); matrix::operator(size_t i, size_t j);`
ベクトルと行列の代入	`vector::operator = (const vector_expression &); vector::operator += (const vector_expression &); vector::operator -= (const vector_expression &); vector::operator = (const scalar_expression &); matrix::operator = (const matrix_expression &); matrix::operator += (const matrix_expression &); matrix::operator -= (const matrix_expression &); matrix::operator = (const scalar_expression &);`
ベクトルと行列に対する単項演算	`vector_expression operator - (const vector_expression &); matrix_expression operator - (const matrix_expression &);`
ベクトルと行列に対する2項演算	`vector_expression operator + (const vector_expression &, const vector_expression &); vector_expression operator - (const vector_expression &, const vector_expression &); matrix_expression operator + (const matrix_expression &, const matrix_expression &); matrix_expression operator - (const matrix_expression &, const matrix_expression &);`
ベクトルと行列の、スカラーとの乗算	`vector_expression operator * (const scalar_expression &, const vector_expression &); vector_expression operator * (const vector_expression &, const scalar_expression &); matrix_expression operator * (const scalar_expression &, const matrix_expression &); matrix_expression operator * (const matrix_expression &, const scalar_expression &);`

次の通り、他の原初的なものに対しては演算子オーバーロードを使わないことにした。

説明	関数
ベクトルが左項のベクトル・行列の乗算	`vector_expression prod<vector_type> (const matrix_expression &, const vector_expression &); vector_expression prod (const matrix_expression &, const vector_expression &);`
ベクトルが右項の行列・ベクトルの乗算	`vector_expression prod<vector_type> (const vector_expression &, const matrix_expression &); vector_expression prod (const vector_expression &, const matrix_expression &);`
行列同士の乗算	`matrix_expression prod<matrix_type> (const matrix_expression &, const matrix_expression &); matrix_expression prod (const matrix_expression &, const matrix_expression &);`
ベクトルの内積	`scalar_expression inner_prod (const vector_expression &, const vector_expression &);`
ベクトル同士の乗算(返り値は行列)	`matrix_expression outer_prod (const vector_expression &, const vector_expression &);`
行列の転置	`matrix_expression trans (const matrix_expression &);`

Efficiency

数値計算の効率という目的を達成するために、 C++ で抽象化を定式化する上でのふたつの困難を克服しなければならない。つまり、一時オブジェクトと仮想関数呼び出しである。式テンプレートはこれらの問題を解決するが、コンパイル時間を遅くする傾向がある。

一時オブジェクトの排除

ベクトルと行列に対する抽象的な式は通常、単項及び2項演算の組を構成する。そのような式を評価する伝統的な方法は、まず合成関数の全ての構成要素を評価して一時オブジェクトにして、続いて合成された結果を評価して、別の一時オブジェクトにするというものである。この方法は、時間効率という点では小さな、空間効率という点では大きなベクトルと行列に対して、高くつく。この問題を解決するアプローチは、最近の関数志向プログラミング言語から知られているように、遅延評価を使うことである。このアプローチの原理は、複雑な式の要素を賢く評価し、対象に直接代入するというものである。

二つの興味深く、危険な事実がある。

まず、ベクトルと行列に対する、要素の賢い評価を使うことの深刻な副作用を被るかもしれない。行列とベクトルの積 x = A x を考えよう。 A₁x の評価と、 x₁ への代入は、右側を変化させるので、 A₂x の評価は誤った結果となる。この問題に対する我々の解決法は、代入の右辺を評価して一時オブジェクトにいれ、それからこの一時オブジェクトを左辺に代入するというものである。さらなる最適化を可能にするために、全ての代入演算子に対応するメンバ関数を提供している。このメンバ関数を使うことで、プログラマは代入の左辺と右辺が独立であり、要素の賢い評価と、対象への直接の代入が安全であることを確かなものとできる。

次に、計算量がひどく悪くなってしまう状況がある。連鎖した行列とベクトルの積 A (B x) を考えてみよう。 A (B x) の従来の評価は二次式である。 B x_i の遅延評価は線形である。 B x_i の全ての要素は、その大きさに依存した線形時間を必要とするので、連鎖した行列とベクトルの積 A (B x) の完全な遅延評価は三次式である。そのような場合には、式の中に一時オブジェクトを再度取り入れる必要がある。

仮想関数呼び出しの排除

式の遅延評価は通常、項のクラス階層の定義となる。これは結果として、ベクトルと行列の単一の要素にアクセスするために動的多相を利用することになるが、動的多相は時間に関して高くつくことが知られている。 2,3年前、David Vandervoorde と Todd Veldhuizen はそれぞれ、一般的には式テンプレートと呼ばれる解決法を発見した。式テンプレートは遅延評価を含み、動的多相を静的、つまりコンパイル時多相に置き換える。式テンプレートは有名な Barton-Nackmann のトリックに強く依存していて、 '奇妙に定義された再帰的テンプレート'という言葉が Jim Coplien によって作られた。

式テンプレートは我々の実装の基礎となっている。

これもよく知られている事実だが、式テンプレートは現在使えるコンパイラに対する挑戦である。 Barton-Nackman のトリックを使うことで結果的に、必要な式テンプレートの量を著しく減らすことが出来た。

また、我々は従来の実装(つまり、式テンプレートを使わない)もサポートしている。これは、開発サイクルをサポートするために、ベクトルと行列の演算についての広範囲の領域限界と型チェックを備えている。デバッグモードからリリースモードへの切り替えは、 <cassert> の NDEBUG プリプロセッサシンボルによって行われている。

Functionality

線形代数をサポートする全ての C++ ライブラリは、長期間存在してきた Fortran のパッケージ BLAS と比較されるだろう。ここでは、 BLAS の呼び出しが我々のクラスにどうマッピングされているかを説明する。

Blas Level 1

BLAS Call	Mapped Library Expression	Mathematical Description	Comment
`_asum`	`norm_1 (x)`	sum \|x_i\|	ベクトルのノルムの合計を計算する。
`_nrm2`	`norm_2 (x)`	sqrt (sum \|x_i\|² )	ベクトルのユークリッドノルムを計算する。
`i_amax`	`norm_inf (x) norm_inf_index (x)`	max \|x_i\|	ベクトルの最大ノルムを計算する。 BLAS はこの値を持つ最初の要素のインデックスを計算している。
`_dot _dotu _dotc`	`inner_prod (x, y)`or`inner_prod (conj (x), y)`	x^T y or x^H y	二つのベクトルの内積を計算する。 BLAS はいくつかのループの展開を実装している。
`dsdot sdsdot`	`a + prec_inner_prod (x, y)`	a + x^T y	double の精度で内積を計算する。
`_copy`	`x = y y.assign (x)`	x <- y	ベクトルを別のベクトルにコピーする。 BLAS はいくつかのループの展開を実装している。
`_swap`	`swap (x, y)`	x <-> y	二つのベクトルを入れ替える。 BLAS はいくつかのループの展開を実装している。
`_scal csscal zdscal`	`x *= a`	x <- a x	ベクトルをスカラー倍する。 BLAS はいくつかのループの展開を実装している。
`_axpy`	`y += a * x`	y <- a x + y	スカラー倍されたベクトルを加える。 BLAS はいくつかのループの展開を実装している。
`_rot _rotm csrot zdrot`	`t.assign (a * x + b * y), y.assign (- b * x + a * y), x.assign (t)`	(x, y) <- (a x + b y, -b x + a y)	平面回転を適用する。
`_rotg _rotmg`		(a, b) <- (? a / sqrt (a² + b²), ? b / sqrt (a² + b²)) or (1, 0) <- (0, 0)	平面回転を構築する。

Blas Level 2

BLAS Call	Mapped Library Expression	Mathematical Description	Comment
`_t_mv`	`x = prod (A, x)` or`x = prod (trans (A), x)` or`x = prod (herm (A), x)`	x <- A x or x <- A^T x or x <- A^H x	行列とベクトルの積を計算する。
`_t_sv`	`y = solve (A, x, tag)` or `inplace_solve (A, x, tag)` or `y = solve (trans (A), x, tag)` or`inplace_solve (trans (A), x, tag)`or`y = solve (herm (A), x, tag)`or`inplace_solve (herm (A), x, tag)`	y <- A^-1 x or x <- A^-1 x or y <- A^T^{^-1} x or x <- A^T^{^-1} x or y <- A^H^{^-1} x or x <- A^H^{^-1} x	三角行列の形式で連立一次方程式を解く。つまり、 A は三角行列である。
`_g_mv _s_mv _h_mv`	`y = a * prod (A, x) + b * y` or`y = a * prod (trans (A), x) + b * y`or`y = a * prod (herm (A), x) + b * y`	y <- a A x + b y or y <- a A^T x + b y y <- a A^H x + b y	行列とベクトルの積のスカラー倍を加える。
`_g_r _g_ru _g_rc`	`A += a * outer_prod (x, y)` or`A += a * outer_prod (x, conj (y))`	A <- a x y^T + A or A <- a x y^H + A	階数 1 の更新を実行する。
`_s_r _h_r`	`A += a * outer_prod (x, x)` or`A += a * outer_prod (x, conj (x))`	A <- a x x^T + A or A <- a x x^H + A	対称行列またはエルミート行列の階数 1 の更新を実行する。
`_s_r2 _h_r2`	`A += a * outer_prod (x, y) + a * outer_prod (y, x))`or`A += a * outer_prod (x, conj (y)) + conj (a) * outer_prod (y, conj (x)))`	A <- a x y^T + a y x^T + A or A <- a x y^H + a^- y x^H + A	対称行列またはエルミート行列の階数 2 の更新を実行する。

Blas Level 3

BLAS Call	Mapped Library Expression	Mathematical Description	Comment
`_t_mm`	`B = a * prod (A, B)` or `B = a * prod (trans (A), B)` or `B = a * prod (A, trans (B))` or `B = a * prod (trans (A), trans (B))` or `B = a * prod (herm (A), B)` or `B = a * prod (A, herm (B))` or `B = a * prod (herm (A), trans (B))` or `B = a * prod (trans (A), herm (B))` or `B = a * prod (herm (A), herm (B))`	B <- a op (A) op (B) with op (X) = X or op (X) = X^T or op (X) = X^H	二つの行列の積のスカラー倍を計算する。
`_t_sm`	`C = solve (A, B, tag)` or `inplace_solve (A, B, tag)` or `C = solve (trans (A), B, tag)` or`inplace_solve (trans (A), B, tag)`or`C = solve (herm (A), B, tag)` or`inplace_solve (herm (A), B, tag)`	C <- A^-1 B or B <- A^-1 B or C <- A^T^{^-1} B or B <- A^-1 B or C <- A^H^{^-1} B or B <- A^H^{^-1} B	三角行列の係数を持つ連立一次方程式を解く。つまり、 A は三角行列である。
`_g_mm _s_mm _h_mm`	`C = a * prod (A, B) + b * C` or `C = a * prod (trans (A), B) + b * C` or `C = a * prod (A, trans (B)) + b * C` or `C = a * prod (trans (A), trans (B)) + b * C` or `C = a * prod (herm (A), B) + b * C` or `C = a * prod (A, herm (B)) + b * C` or `C = a * prod (herm (A), trans (B)) + b * C` or `C = a * prod (trans (A), herm (B)) + b * C` or `C = a * prod (herm (A), herm (B)) + b * C`	C <- a op (A) op (B) + b C with op (X) = X or op (X) = X^T or op (X) = X^H	二つの行列の積のスカラー倍を加える。
`_s_rk _h_rk`	`B = a * prod (A, trans (A)) + b * B` or `B = a * prod (trans (A), A) + b * B` or `B = a * prod (A, herm (A)) + b * B` or `B = a * prod (herm (A), A) + b * B`	B <- a A A^T + b B or B <- a A^T A + b B or B <- a A A^H + b B or B <- a A^H A + b B	対称行列またはエルミート行列の階数 k の更新を実行する。
`_s_r2k _h_r2k`	`C = a * prod (A, trans (B)) + a * prod (B, trans (A)) + b * C`or `C = a * prod (trans (A), B) + a * prod (trans (B), A) + b * C`or `C = a * prod (A, herm (B)) + conj (a) * prod (B, herm (A)) + b * C`or `C = a * prod (herm (A), B) + conj (a) * prod (herm (B), A) + b * C`	C <- a A B^T + a B A^T + b C or C <- a A^T B + a B^T A + b C or C <- a A B^H + a^- B A^H + b C or C <- a A^H B + a^- B^H A + b C	対称行列またはエルミート行列の階数 2 k の更新を実行する。

Storage Layouts

ライブラリは伝統的な密ベクトルと密行列、圧縮されたベクトルと行列、基本的な疎ベクトルと疎行列といった記憶域レイアウトをサポートしている。ベクトルと行列の最も一般的な構築の記述は以下の通りである。

Construction	Comment
`vector<T, std::vector<T> > v (size)`	密ベクトルの構築。記憶域は標準のベクタ(std::vector)によって提供される。記憶域レイアウトは通常、 BLAS に対応している。
`vector<T, unbounded_array<T> > v (size)`	密ベクトルの構築。記憶域はヒープに作られる配列によって提供される。記憶域レイアウトは通常、 BLAS に対応している。
`vector<T, bounded_array<T, N> > v (size)`	Constructs a dense vector, storage is provided by a stack-based array. The storage layout usually is BLAS compliant. 記憶域レイアウトは通常、 BLAS に対応している。
`unit_vector<T> v (size, index)`	`index` 次の典型的な単位ベクトルの構築。
`zero_vector<T> v (size)`	ゼロベクトルの構築。
`sparse_vector<T, std::map<std::size_t, T> > v (size, non_zeros)`	疎ベクトルの構築。記憶域は標準のマップ(std::map)によって提供される。
`sparse_vector<T, map_array<std::size_t, T> > v (size, non_zeros)`	疎ベクトルの構築。記憶域はマップ配列によって提供される。
`compressed_vector<T> v (size, non_zeros)`	圧縮行列の構築。記憶域レイアウトは通常、 BLAS に対応している。
`coordinate_vector<T> v (size, non_zeros)`	座標行列の構築。記憶域レイアウトは通常、 BLAS に対応している。
`vector_range<V> vr (v, range)`	密ベクトル、圧縮されたベクトル、疎ベクトルの部分ベクトルを構築する。
`vector_slice<V> vs (v, slice)`	密ベクトル、圧縮されたベクトル、疎ベクトルの部分ベクトルを構築する。このクラスは通常、 BLAS の vector slices をエミュレートするために使われる。
`matrix<T, row_major, std::vector<T> > m (size1, size2)`	密行列の構築。方向は行優先である。記憶域は標準のベクタによって提供される。
`matrix<T, column_major, std::vector<T> > m (size1, size2)`	密行列の構築。方向は列優先である。記憶域は標準のベクタによって提供される。記憶域レイアウトは通常、 BLAS に対応している。
`matrix<T, row_major, unbounded_array<T> > m (size1, size2)`	密行列の構築。方向は行優先である。記憶域はヒープに作られる配列によって提供される。
`matrix<T, column_major, unbounded_array<T> > m (size1, size2)`	密行列の構築。方向は列優先である。記憶域はヒープに作られる配列によって提供される。記憶域レイアウトは通常、 BLAS に対応している。
`matrix<T, row_major, bounded_array<T, N1 * N2> > m (size1, size2)`	密行列の構築。方向は行優先である。記憶域はスタックに作られる配列によって提供される。
`matrix<T, column_major, bounded_array<T, N1 * N2> > m (size1, size2)`	密行列の構築。方向は列優先である。記憶域はスタックに作られる配列によって提供される。記憶域レイアウトは通常、 BLAS に対応している。
`identity_matrix<T> m (size)`	単位行列の構築。
`zero_matrix<T> m (size1, size2)`	ゼロ行列の構築。
`triangular_matrix<T, row_major, F, A> m (size)`	圧縮された三角行列の構築。方向は行優先である。
`triangular_matrix<T, column_major, F, A> m (size)`	圧縮された三角行列の構築。方向は列優先である。記憶域レイアウトは通常、 BLAS に対応している。
`banded_matrix<T, row_major, A> m (size1, size2, lower, upper)`	圧縮された帯行列の構築。方向は行優先である。
`banded_matrix<T, column_major, A> m (size1, size2, lower, upper)`	圧縮された帯行列の構築。方向は列優先である。記憶域レイアウトは通常、 BLAS に対応している。
`symmetric_matrix<T, row_major, F, A> m (size)`	圧縮された対称行列の構築。方向は行優先である。
`symmetric_matrix<T, column_major, F, A> m (size)`	圧縮された対称行列の構築。方向は列優先である。記憶域レイアウトは通常、 BLAS に対応している。
`hermitian_matrix<T, row_major, F, A> m (size)`	圧縮されたエルミート行列の構築。方向は行優先である。
`hermitian_matrix<T, column_major, F, A> m (size)`	圧縮されたエルミート行列の構築。方向は列優先である。記憶域レイアウトは通常、 BLAS に対応している。
`sparse_matrix<T, row_major, std::map<std::size_t, T> > m (size1, size2, non_zeros)`	疎行列の構築。方向は行優先である。記憶域は標準のマップによって提供される。
`sparse_matrix<T, column_major, std::map<std::size_t, T> > m (size1, size2, non_zeros)`	疎行列の構築。方向は列優先である。記憶域は標準のマップによって提供される。
`sparse_matrix<T, row_major, map_array<std::size_t, T> > m (size1, size2, non_zeros)`	疎行列の構築。方向は行優先である。記憶域はマップ配列によって提供される。
`sparse_matrix<T, column_major, map_array<std::size_t, T> > m (size1, size2, non_zeros)`	疎行列の構築。方向は列優先である。記憶域はマップ配列によって提供される。
`compressed_matrix<T, row_major> m (size1, size2, non_zeros)`	圧縮行列の構築。方向は行優先である。記憶域レイアウトは通常、 BLAS に対応している。
`compressed_matrix<T, column_major> m (size1, size2, non_zeros)`	圧縮行列の構築。方向は列優先である。記憶域レイアウトは通常、 BLAS に対応している。
`coordinate_matrix<T, row_major> m (size1, size2, non_zeros)`	座標行列の構築。方向は行優先である。記憶域レイアウトは通常、 BLAS に対応している。
`coordinate_matrix<T, column_major> m (size1, size2, non_zeros)`	座標行列の構築。方向は列優先である。記憶域レイアウトは通常、 BLAS に対応している。
`matrix_row<M> mr (m, i)`	密行列、圧縮された行列、疎行列の、i 行目を含む部分ベクトルの構築。
`matrix_column<M> mc (m, j)`	密行列、圧縮された行列、疎行列の、j 列目を含む部分ベクトルの構築。
`matrix_range<M> mr (m, range1, range2)`	密行列、圧縮された行列、疎行列の部分ベクトルの構築。このクラスは通常、 BLAS のleading dimensions をエミュレートするために使われる。
`matrix_slice<M> ms (m, slice1, slice2)`	密行列、圧縮された行列、疎行列の部分ベクトルの構築。
`triangular_adaptor<M, F> ta (m)`	密行列、圧縮された行列、疎行列からの三角行列の構築。このクラスは通常、対応する BLAS の行列型を生成するために使われる。
`banded_adaptor<M> ba (m, lower, upper)`	密行列、圧縮された行列、疎行列からの帯行列の構築。このクラスは通常、対応する BLAS の行列型を生成するために使われる。
`symmetric_adaptor<M> sa (m)`	密行列、圧縮された行列、疎行列からの対称行列の構築。このクラスは通常、対応する BLAS の行列型を生成するために使われる。
`hermitian_adaptor<M> ha (m)`	密行列、圧縮された行列、疎行列からのエルミート行列の構築。このクラスは通常、対応する BLAS の行列型を生成するために使われる。

Compatibility

互換性という理由で、ベクトルと行列に対して配列風のインデックス付けを提供している。これはしかし、一時プロキシオブジェクトが必要なために、行列に対して高くつくこともある。

最も広く使われている C++ コンパイラをサポートするために、我々の設計と実装は基本的に、テンプレートの部分特殊化には依存していない。

ライブラリは operator new と operator delete による、標準準拠のメモリ割り当てを仮定している。そのため、 MSVC 6.0 で実行するように意図されたプログラムはメモリの状態を検出するために、 _set_new_handler によって、 std::bad_alloc 例外を投げる正しい new ハンドラを設定するべきである。

MSVC 6.0 で気が狂うほどの効率を得るためには、ヘッダファイル config.hppの中で、BOOST_UBLAS_INLINE のプリプロセッサ定義を __forceinline に変更すべきである。しかし我々は、この最適化が脆弱なものではないかと疑っている。

Reference

Benchmark Results

The following tables contain results of one of our benchmarks. This benchmark compares a native C implementation ('C array') and some library based implementations. The safe variants based on the library assume aliasing, the fast variants do not use temporaries and are functionally equivalent to the native C implementation. Besides the generic vector and matrix classes the benchmark utilizes special classes c_vector and c_matrix, which are intended to avoid every overhead through genericity.

The benchmark program was compiled with MSVC 6.0 and run on an Intel Pentium II with 333 Mhz. First we comment the results for double vectors and matrices of dimension 3 and 3 x 3, respectively.

Operation	Implementation	Elapsed [s]	MFLOP/s	Comment
inner_prod	C array	0.1	47.6837	No significant abstraction penalty
	c_vector	0.06	79.4729
	vector<unbounded_array>	0.11	43.3488
	vector<std::vector>	0.11	43.3488
vector + vector	C array	0.05	114.441	Abstraction penalty: factor 2
	c_vector safe	0.22	26.0093
	c_vector fast	0.11	52.0186
	vector<unbounded_array> safe	1.05	5.44957
	vector<unbounded_array> fast	0.16	35.7628
	vector<std::vector> safe	1.16	4.9328
	vector<std::vector> fast	0.16	35.7628
outer_prod	C array	0.06	85.8307	Abstraction penalty: factor 2
	c_matrix, c_vector safe	0.22	23.4084
	c_matrix, c_vector fast	0.11	46.8167
	matrix<unbounded_array>, vector<unbounded_array> safe	0.38	13.5522
	matrix<unbounded_array>, vector<unbounded_array> fast	0.16	32.1865
	matrix<std::vector>, vector<std::vector> safe	0.5	10.2997
	matrix<std::vector>, vector<std::vector> fast	0.11	46.8167
prod (matrix, vector)	C array	0.06	71.5256	No significant abstraction penalty
	c_matrix, c_vector safe	0.11	39.0139
	c_matrix, c_vector fast	0.11	39.0139
	matrix<unbounded_array>, vector<unbounded_array> safe	0.33	13.0046
	matrix<unbounded_array>, vector<unbounded_array> fast	0.11	39.0139
	matrix<std::vector>, vector<std::vector> safe	0.38	11.2935
	matrix<std::vector>, vector<std::vector> fast	0.05	85.8307
matrix + matrix	C array	0.11	46.8167	No significant abstraction penalty
	c_matrix safe	0.17	30.2932
	c_matrix fast	0.11	46.8167
	matrix<unbounded_array> safe	0.44	11.7042
	matrix<unbounded_array> fast	0.16	32.1865
	matrix<std::vector> safe	0.6	8.58307
	matrix<std::vector> fast	0.17	30.2932
prod (matrix, matrix)	C array	0.11	39.0139	No significant abstraction penalty
	c_matrix safe	0.11	39.0139
	c_matrix fast	0.11	39.0139
	matrix<unbounded_array> safe	0.22	19.507
	matrix<unbounded_array> fast	0.11	39.0139
	matrix<std::vector> safe	0.27	15.8946
	matrix<std::vector> fast	0.11	39.0139

We notice a twofold performance loss for small vectors and matrices: first the general abstraction penalty for using classes, and then a small loss when using the generic vector and matrix classes. The difference w.r.t. alias assumptions is also significant.

Next we comment the results for double vectors and matrices of dimension 100 and 100 x 100, respectively.

Operation	Implementation	Elapsed [s]	MFLOP/s	Comment
inner_prod	C array	0.05	113.869	No significant abstraction penalty
	c_vector	0.06	94.8906
	vector<unbounded_array>	0.05	113.869
	vector<std::vector>	0.06	94.8906
vector + vector	C array	0.05	114.441	No significant abstraction penalty
	c_vector safe	0.11	52.0186
	c_vector fast	0.11	52.0186
	vector<unbounded_array> safe	0.11	52.0186
	vector<unbounded_array> fast	0.06	95.3674
	vector<std::vector> safe	0.17	33.6591
	vector<std::vector> fast	0.11	52.0186
outer_prod	C array	0.05	114.441	No significant abstraction penalty
	c_matrix, c_vector safe	0.28	20.4359
	c_matrix, c_vector fast	0.11	52.0186
	matrix<unbounded_array>, vector<unbounded_array> safe	0.27	21.1928
	matrix<unbounded_array>, vector<unbounded_array> fast	0.06	95.3674
	matrix<std::vector>, vector<std::vector> safe	0.28	20.4359
	matrix<std::vector>, vector<std::vector> fast	0.11	52.0186
prod (matrix, vector)	C array	0.11	51.7585	No significant abstraction penalty
	c_matrix, c_vector safe	0.11	51.7585
	c_matrix, c_vector fast	0.05	113.869
	matrix<unbounded_array>, vector<unbounded_array> safe	0.11	51.7585
	matrix<unbounded_array>, vector<unbounded_array> fast	0.06	94.8906
	matrix<std::vector>, vector<std::vector> safe	0.1	56.9344
	matrix<std::vector>, vector<std::vector> fast	0.06	94.8906
matrix + matrix	C array	0.22	26.0093	No significant abstraction penalty
	c_matrix safe	0.49	11.6776
	c_matrix fast	0.22	26.0093
	matrix<unbounded_array> safe	0.39	14.6719
	matrix<unbounded_array> fast	0.22	26.0093
	matrix<std::vector> safe	0.44	13.0046
	matrix<std::vector> fast	0.27	21.1928
prod (matrix, matrix)	C array	0.06	94.8906	No significant abstraction penalty
	c_matrix safe	0.06	94.8906
	c_matrix fast	0.05	113.869
	matrix<unbounded_array> safe	0.11	51.7585
	matrix<unbounded_array> fast	0.17	33.4908
	matrix<std::vector> safe	0.11	51.7585
	matrix<std::vector> fast	0.16	35.584

For larger vectors and matrices the general abstraction penalty for using classes seems to decrease, the small loss when using generic vector and matrix classes seems to remain. The difference w.r.t. alias assumptions remains visible, too.

Copyright (©) 2000-2002 Joerg Walter, Mathias Koch
Permission to copy, use, modify, sell and distribute this document is granted provided this copyright notice appears in all copies. This document is provided ``as is'' without express or implied warranty, and with no claim as to its suitability for any purpose.

Last revised: 1/15/2003