良い出発点は、Robert A. van de Geijn と Enrique S. Quintana-Ortí による素晴らしい本The Science of Programming Matrix Computationsです。彼らは無料のダウンロード版を提供しています。
BLAS は 3 つのレベルに分かれています。
レベル 1 は、ベクトルのみを操作する一連の線形代数関数を定義します。これらの関数は、ベクトル化 (たとえば、SSE の使用) の恩恵を受けます。
レベル 2 の関数は、行列とベクトルの積などの行列とベクトルの演算です。これらの関数は、レベル 1 関数の観点から実装できます。ただし、共有メモリを備えたマルチプロセッサ アーキテクチャを利用する専用の実装を提供できれば、この関数のパフォーマンスを向上させることができます。
レベル 3 の関数は、行列-行列の積のような操作です。繰り返しますが、Level2 関数に関してそれらを実装できます。しかし、Level3 関数は O(N^2) データに対して O(N^3) 操作を実行します。したがって、プラットフォームにキャッシュ階層がある場合、キャッシュ最適化/キャッシュ フレンドリーな専用の実装を提供すると、パフォーマンスを向上させることができます。これは本の中でうまく説明されています。Level3 機能の主な向上は、キャッシュの最適化によるものです。このブーストは、並列処理やその他のハードウェアの最適化による 2 番目のブーストを大幅に上回ります。
ところで、高性能 BLAS 実装のほとんど (またはすべて) は、Fortran では実装されていません。ATLAS は C で実装されています。GotoBLAS/OpenBLAS は C で実装され、パフォーマンスに重要な部分はアセンブラーで実装されています。BLAS の参照実装のみが Fortran で実装されています。ただし、これらすべての BLAS 実装は、LAPACK とリンクできるように Fortran インターフェイスを提供します (LAPACK は BLAS からすべてのパフォーマンスを得ます)。
最適化されたコンパイラは、この点でマイナーな役割を果たします (GotoBLAS/OpenBLAS の場合、コンパイラはまったく問題になりません)。
IMHO no BLAS の実装では、Coppersmith–Winograd アルゴリズムや Strassen アルゴリズムなどのアルゴリズムを使用しています。考えられる理由は次のとおりです。
- おそらく、これらのアルゴリズムのキャッシュ最適化された実装を提供することは不可能です (つまり、勝つよりも負けるほうが多い)
- これらのアルゴリズムは数値的に安定していません。BLAS は LAPACK の計算カーネルであるため、これは不可能です。
- これらのアルゴリズムは紙の上ではかなりの時間の複雑さを持っていますが、Big O 表記は大きな定数を隠しているため、非常に大きな行列に対してのみ実行可能になり始めます。
編集/更新:
このトピックの新しい画期的な論文はBLIS 論文です。それらは非常によく書かれています。私の講義「ハイ パフォーマンス コンピューティングのためのソフトウェアの基礎」では、彼らの論文に従ってマトリックス-マトリックス製品を実装しました。実際、私はマトリックス-マトリックス製品のいくつかの変形を実装しました。最も単純な亜種はすべてプレーン C で記述されており、コードは 450 行未満です。他のすべてのバリアントは単にループを最適化するだけです
for (l=0; l<MR*NR; ++l) {
AB[l] = 0;
}
for (l=0; l<kc; ++l) {
for (j=0; j<NR; ++j) {
for (i=0; i<MR; ++i) {
AB[i+j*MR] += A[i]*B[j];
}
}
A += MR;
B += NR;
}
マトリックス - マトリックス製品の全体的なパフォーマンスは、これらのループのみに依存します。時間の約 99.9% がここで費やされます。他のバリアントでは、組み込み関数とアセンブラー コードを使用してパフォーマンスを向上させました。ここで、すべてのバリアントを通過するチュートリアルを見ることができます。
ulmBLAS: GEMM (Matrix-Matrix Product) のチュートリアル
BLIS の論文と合わせて、インテル® MKL のようなライブラリーがどのようにそのようなパフォーマンスを達成できるかを理解するのはかなり簡単になります。そして、行または列の主要なストレージを使用するかどうかが問題にならないのはなぜですか!
最終的なベンチマークは次のとおりです (プロジェクトを ulmBLAS と呼びます)。
ulmBLAS、BLIS、MKL、openBLAS、Eigen のベンチマーク
別の編集/更新:
BLAS が線形方程式系を解くような数値線形代数の問題にどのように使用されるかについてのチュートリアルも書きました。
高性能 LU 因数分解
(この LU 因数分解は、たとえば、線形方程式系を解くために Matlab で使用されます。)
時間を見つけたいと思っています。
わかりました、どうぞ:キャッシュ最適化並列 LU 因数分解のコーディング
PS: uBLAS のパフォーマンスを改善するための実験もいくつか行いました。実際、uBLAS のパフォーマンスを向上させるのは非常に簡単です (そうです、言葉遊びです :))。
uBLAS での実験。
BLAZEを使用した同様のプロジェクト:
BLAZE の実験。