[解決済み] BLASはどうしてあんなに高性能なのか？

質問

好奇心から、私は自分自身の行列乗算関数とBLAS実装のベンチマークを取ることにしました... 私は控えめに言っても、その結果に驚きました。

カスタム実装、10回の試行 1000x1000 の行列の乗算を 10 回試行しました。

Took: 15.76542 seconds.

BLASの実装、10回の試行で 1000x1000の行列の乗算を10回試行。

Took: 1.32432 seconds.

これは単精度浮動小数点数を使用しています。

私の実装です。

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

2つ質問があります。

1. 行列と行列の乗算、例えば nxm * mxn が n*n*m 回の乗算を必要とするとすると、上記のケースでは 1000^3 または 1e9 回の演算が必要です。私の2.6Ghzのプロセッサで、BLASが10*1e9の演算を1.32秒で行うことは可能なのでしょうか？乗算が1つの演算で、他に何もしていないとしても、～4秒かかるはずです。
2. なぜ私の実装はこれほど遅いのですか?

どのように解決するのですか？

良い出発点は、偉大な本である 行列計算のプログラミングの科学 Robert A. van de Geijn と Enrique S. Quintana-Ortí によるものです。Quintana-Ortí です。 無料ダウンロード版もあります。

BLASは3つのレベルに分かれています。

• レベル1では、ベクトルに対してのみ演算する線形代数関数のセットを定義しています。 これらの関数はベクトル化（例えばSSEの使用）から利益を得ます。

• レベル2関数は、行列-ベクトル演算、例えば、いくつかの行列-ベクトル積です。 これらの関数は、レベル1関数で実装することができます。 しかし、共有メモリを持つ何らかのマルチプロセッサアーキテクチャを利用した専用の実装を提供できれば、この関数の性能を高めることができます。

• レベル3関数は、行列と行列の積のような演算です。 ここでもレベル2関数で実装することができます。 しかし、レベル3関数はO(N^2)個のデータに対してO(N^3)個の演算を実行します。 ですから、もしあなたのプラットフォームがキャッシュ階層を持つなら、以下のような専用の実装を提供すれば、性能を向上させることができます。 キャッシュ最適化/キャッシュフレンドリー . これについては、この本でうまく説明されています。 Level3 関数の主なブーストはキャッシュの最適化から生まれます。 このブーストは、並列処理やその他のハードウェア最適化による2番目のブーストを大幅に上回ります。

ところで、高性能な BLAS 実装のほとんど (あるいはすべて) は Fortran で実装されているわけではありません。 ATLAS は C で実装され、GotoBLAS/OpenBLAS は C で、その性能的に重要な部分は Assembler で実装されています。 BLASのリファレンス実装のみがFortranで実装されています。 しかし、これらのBLASの実装はすべて、LAPACKとリンクできるようにFortranインターフェースを提供しています（LAPACKはBLASからすべての性能を獲得しています）。

最適化されたコンパイラーはこの点で小さな役割を果たします（そして GotoBLAS/OpenBLAS ではコンパイラーはまったく重要ではありません）。

IMHOでは、Coppersmith-WinogradアルゴリズムやStrassenアルゴリズムのようなアルゴリズムを使用するBLAS実装はありません。 考えられる理由は

• これらのアルゴリズムのキャッシュ最適化された実装を提供することが不可能なため (つまり、勝つことよりも失うことの方が多い)。
• これらのアルゴリズムは数値的に安定していません。 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 (行列-行列積)のチュートリアル

BLISの論文と合わせると、Intel MKLのようなライブラリがなぜこのようなパフォーマンスを得られるのか、かなり簡単に理解できるようになります。そして、行または列の主要なストレージを使用するかどうかは、なぜ重要ではないのでしょうか!

最終的なベンチマークはこちらです (私たちはこのプロジェクトを ulmBLAS と呼んでいます)。

ulmBLAS、BLIS、MKL、openBLAS、Eigenのベンチマークです。

もう一つの編集・更新。

BLASが連立一次方程式の解のような数値線形代数問題にどのように使用されるかについてのチュートリアルも書きました。

高性能LU因子分解

(このLU因数分解は、例えばMatlabで連立一次方程式の解法に使われます)。

のようなLU因数分解の高度にスケーラブルな並列実装を実現する方法を説明し実証するために、このチュートリアルを拡張する時間があればと思います。 PLASMA .

よし、これでいいんだ。 キャッシュ最適化並列LU因数分解のコード化

追記: uBLASの性能向上に関する実験も行いました。 実はuBLASの性能を上げるのはとても簡単です（言葉遊びです :)。

uBLASの実験 .

ここでは、同様のプロジェクトで BLAZE :

BLAZEに関する実験 .