2つの3x3行列の積を取りA*B=Cます。単純に、これには標準アルゴリズムを使用した27回の乗算が必要です。賢い人なら、23回の掛け算だけでこれを行うことができます。これは1973年にLadermanによって発見された結果です。この手法では、中間ステップを保存し、それらを適切な方法で組み合わせる必要があります。
次に、言語と型を修正しましょう。たとえば、の要素を持つC++ですdouble。Ladermanアルゴリズムが単純な二重ループに対してハードコーディングされている場合、最新のコンパイラーのパフォーマンスがアルゴリズムの違いを打ち消すことを期待できますか？
この質問に関する注記：これはプログラミングサイトであり、タイムクリティカルな内部ループのベストプラクティスのコンテキストで質問が行われます。時期尚早の最適化ではありません。実装のヒントはコメントとして大歓迎です。
重要なのは、プラットフォームで設定された命令をマスターすることです。プラットフォームによって異なります。いくつかの手法があり、可能な限り最大のパフォーマンスが必要になる傾向がある場合、コンパイラにはプロファイリングツールが付属し、その一部には最適化のヒントが組み込まれています。最も細かい操作については、アセンブラの出力を調べて、改善があるかどうかを確認してください。そのレベルでも。
同時命令複数のデータコマンドは、複数のオペランドに対して同じ操作を並行して実行します。あなたが取ることができるように
double a,b,c,d;
double w = d + a;
double x = a + b;
double y = b + c;
double z = c + d;
と置き換えます
double256 dabc = pack256(d, a, b, c);
double256 abcd = pack256(a, b, c, d);
double256 wxyz = dabc + abcd;
そのため、値がレジスタにロードされると、256ビット幅のレジスタを備えた架空のプラットフォームの単一の256ビット幅のレジスタにロードされます。
浮動小数点は重要な考慮事項です。一部のDSPは、整数での動作が大幅に高速になる可能性があります。GPUは浮動小数点で優れている傾向がありますが、単精度で2倍高速なものもあります。この問題の3x3のケースは、単一のCUDAスレッドに収まる可能性があるため、これらの計算のうち256を同時にストリーミングできます。
プラットフォームを選択し、ドキュメントを読み、いくつかの異なるメソッドを実装して、それらのプロファイルを作成します。
私は自分でタイミングテストを実行しましたが、その結果には驚かされました（そのため、最初に質問したのはなぜですか）。要するに、標準のコンパイルでladermanは約225％速くなりますが、-03最適化フラグを使用すると50％遅くなります。フラグを立てるたびにランダムな要素をマトリックスに追加する必要がありました。そうし-O3ないと、コンパイラーは単純な乗算を完全に最適化し、クロック精度の範囲内でゼロの時間を取りました。ladermanアルゴリズムはチェック/ダブルチェックするのが面倒だったので、後世のために以下の完全なコードを投稿します。
仕様：Ubuntu 12.04、Dell Prevision T1600、gcc。タイミングのパーセント差：
g++ [2.22, 2.23, 2.27] g++ -O3 [-0.48, -0.49, -0.48]g++ -funroll-loops -O3 [-0.48, -0.48, -0.47] Ladermanの実装に伴うベンチマークコード：
#include <iostream>
#include <ctime>
#include <cstdio>
#include <cstdlib>
using namespace std;
void simple_mul(const double a[3][3],
const double b[3][3],
double c[3][3]) {
int i,j,m,n;
for(i=0;i<3;i++) {
for(j=0;j<3;j++) {
c[i][j] = 0;
for(m=0;m<3;m++)
c[i][j] += a[i][m]*b[m][j];
}
}
}
void laderman_mul(const double a[3][3],
const double b[3][3],
double c[3][3]) {
double m[24]; // not off by one, just wanted to match the index from the paper
m[1 ]= (a[0][0]+a[0][1]+a[0][2]-a[1][0]-a[1][1]-a[2][1]-a[2][2])*b[1][1];
m[2 ]= (a[0][0]-a[1][0])*(-b[0][1]+b[1][1]);
m[3 ]= a[1][1]*(-b[0][0]+b[0][1]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][2]);
m[4 ]= (-a[0][0]+a[1][0]+a[1][1])*(b[0][0]-b[0][1]+b[1][1]);
m[5 ]= (a[1][0]+a[1][1])*(-b[0][0]+b[0][1]);
m[6 ]= a[0][0]*b[0][0];
m[7 ]= (-a[0][0]+a[2][0]+a[2][1])*(b[0][0]-b[0][2]+b[1][2]);
m[8 ]= (-a[0][0]+a[2][0])*(b[0][2]-b[1][2]);
m[9 ]= (a[2][0]+a[2][1])*(-b[0][0]+b[0][2]);
m[10]= (a[0][0]+a[0][1]+a[0][2]-a[1][1]-a[1][2]-a[2][0]-a[2][1])*b[1][2];
m[11]= a[2][1]*(-b[0][0]+b[0][2]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][1]);
m[12]= (-a[0][2]+a[2][1]+a[2][2])*(b[1][1]+b[2][0]-b[2][1]);
m[13]= (a[0][2]-a[2][2])*(b[1][1]-b[2][1]);
m[14]= a[0][2]*b[2][0];
m[15]= (a[2][1]+a[2][2])*(-b[2][0]+b[2][1]);
m[16]= (-a[0][2]+a[1][1]+a[1][2])*(b[1][2]+b[2][0]-b[2][2]);
m[17]= (a[0][2]-a[1][2])*(b[1][2]-b[2][2]);
m[18]= (a[1][1]+a[1][2])*(-b[2][0]+b[2][2]);
m[19]= a[0][1]*b[1][0];
m[20]= a[1][2]*b[2][1];
m[21]= a[1][0]*b[0][2];
m[22]= a[2][0]*b[0][1];
m[23]= a[2][2]*b[2][2];
c[0][0] = m[6]+m[14]+m[19];
c[0][1] = m[1]+m[4]+m[5]+m[6]+m[12]+m[14]+m[15];
c[0][2] = m[6]+m[7]+m[9]+m[10]+m[14]+m[16]+m[18];
c[1][0] = m[2]+m[3]+m[4]+m[6]+m[14]+m[16]+m[17];
c[1][1] = m[2]+m[4]+m[5]+m[6]+m[20];
c[1][2] = m[14]+m[16]+m[17]+m[18]+m[21];
c[2][0] = m[6]+m[7]+m[8]+m[11]+m[12]+m[13]+m[14];
c[2][1] = m[12]+m[13]+m[14]+m[15]+m[22];
c[2][2] = m[6]+m[7]+m[8]+m[9]+m[23];
}
int main() {
int N = 1000000000;
double A[3][3], C[3][3];
std::clock_t t0,t1;
timespec tm0, tm1;
A[0][0] = 3/5.; A[0][1] = 1/5.; A[0][2] = 2/5.;
A[1][0] = 3/7.; A[1][1] = 1/7.; A[1][2] = 3/7.;
A[2][0] = 1/3.; A[2][1] = 1/3.; A[2][2] = 1/3.;
t0 = std::clock();
for(int i=0;i<N;i++) {
// A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
simple_mul(A,A,C);
}
t1 = std::clock();
double tdiff_simple = (t1-t0)/1000.;
cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl;
cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl;
cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl;
cout << tdiff_simple << endl;
cout << endl;
t0 = std::clock();
for(int i=0;i<N;i++) {
// A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
laderman_mul(A,A,C);
}
t1 = std::clock();
double tdiff_laderman = (t1-t0)/1000.;
cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl;
cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl;
cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl;
cout << tdiff_laderman << endl;
cout << endl;
double speedup = (tdiff_simple-tdiff_laderman)/tdiff_laderman;
cout << "Approximate speedup: " << speedup << endl;
return 0;
}
主なパフォーマンスの懸念はメモリの待ち時間だと思います。Adouble[9]は通常72バイトです。それはすでに取るに足らない量であり、あなたはそれらのうちの3つを使用しています。
質問ではC++について言及しましたが、C＃（.NET 4.5）で3x3行列乗算C = A * Bを実装し、最適化を行った64ビットWindows7マシンでいくつかの基本的なタイミングテストを実行しました。10,000,000回の乗算に約
興味深いことに、ladermanコードは素朴な方法よりも遅かった。私はプロファイラーで調査しませんでしたが、追加の割り当ては、いくつかの追加の乗算よりもコストがかかると思います。
現在のコンパイラは、これらの最適化を実行するのに十分賢いようです。これは良いことです。これが私があなたの興味のために使った素朴なコードです：
public static Matrix3D operator *(Matrix3D a, Matrix3D b)
{
double c11 = a.M11 * b.M11 + a.M12 * b.M21 + a.M13 * b.M31;
double c12 = a.M11 * b.M12 + a.M12 * b.M22 + a.M13 * b.M32;
double c13 = a.M11 * b.M13 + a.M12 * b.M23 + a.M13 * b.M33;
double c21 = a.M21 * b.M11 + a.M22 * b.M21 + a.M23 * b.M31;
double c22 = a.M21 * b.M12 + a.M22 * b.M22 + a.M23 * b.M32;
double c23 = a.M21 * b.M13 + a.M22 * b.M23 + a.M23 * b.M33;
double c31 = a.M31 * b.M11 + a.M32 * b.M21 + a.M33 * b.M31;
double c32 = a.M31 * b.M12 + a.M32 * b.M22 + a.M33 * b.M32;
double c33 = a.M31 * b.M13 + a.M32 * b.M23 + a.M33 * b.M33;
return new Matrix3D(
c11, c12, c13,
c21, c22, c23,
c31, c32, c33);
}
ここで、Matrix3Dは不変の構造体（読み取り専用のdoubleフィールド）です。
トリッキーなことは、コードを測定するのではなく、コンパイラがコードを使って何をしたかではなく、有効なベンチマークを考え出すことです（大量の余分なものを含むデバッガー、または結果が使用されなかったため、実際のコードなしで最適化されます）。私は通常、コンパイラがテスト対象のコードを削除できないように、結果を「タッチ」しようとします（たとえば、行列要素が89038.8989384と等しいかどうかを確認し、等しい場合はスローします）。しかし、結局、コンパイラがこの比較を邪魔にならないようにハックするかどうかさえわかりません:)