BLAS はどのようにしてこのような優れたパフォーマンスを実現しているのでしょうか? 質問する

良い出発点は素晴らしい本です行列計算プログラミングの科学著者：Robert A. van de Geijn、Enrique S. Quintana-Ortí。無料ダウンロード版も提供しています。

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

• レベル 1 は、ベクトルのみを操作する線形代数関数のセットを定義します。これらの関数は、ベクトル化 (たとえば、SSE などの SIMD の使用) の恩恵を受けます。

• レベル 2 関数は、行列ベクトル演算 (行列ベクトル積など) です。これらの関数は、レベル 1 関数に基づいて実装できます。ただし、共有メモリを備えたマルチプロセッサ アーキテクチャを使用する専用の実装を提供できれば、これらの関数のパフォーマンスを向上させることができます。

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

ちなみに、高性能 BLAS 実装のほとんど (またはすべて) は Fortran で実装されていません。ATLAS は C で実装されています。GotoBLAS/OpenBLAS は C で実装され、パフォーマンスが重要な部分はアセンブラで実装されています。BLAS のリファレンス実装のみが Fortran で実装されています。ただし、これらの BLAS 実装はすべて Fortran インターフェイスを提供しているため、LAPACK にリンクできます (LAPACK は BLAS からすべてのパフォーマンスを得ています)。

最適化されたコンパイラは、この点では小さな役割を果たします (GotoBLAS/OpenBLAS の場合、コンパイラはまったく重要ではありません)。

私の意見では、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;
    }