ほとんどの場合、記事全体を読むにはあまりにも面倒ですが、下にスクロールして結論を読むことをお勧めします-興味深い数字を理解できます。
そのため、タスク：2つの大きなdouble行列（3次サイズ）を乗算します。 すべてのアルゴリズムは長方形のものにも適していますが、簡単にするために、正方形の行列を考えます。 アルゴリズムはC ++で作成されましたが、どこでもクラスを使用したことがないため、コードはC互換と見なすことができます（coutのみが使用した可能性があります）。
ここでは、マトリックスとは何か、それらを乗算する方法については説明しません-これを知らない人は、乗算を加速する方法にほとんど興味がありません...
まず、行列がサイズn * nの1次元配列に格納されることを決定します。 要素A（i、j）は[i * n + j]と呼ばれます。 2次元配列よりも高速に1次元に格納するため、 1つの整数乗算と1つのメモリアクセスは、2つのメモリアクセスよりも高速です。
乗算する最初の最も論理的な方法
int multMatrixSqBad（double * a、const int n、double * b、double * c）{
int i、j、k;
for（i = 0; i <n; i ++）{
for（j = 0; j <n; j ++）{
c [i * n + j] = 0;
for（k = 0; k <n; k ++）{
c [i * n + j] + = a [i * n + k] * b [k * n + j];
}
}
}
0を返します。
}
最初に、サイクルでi * nを何回も掛けすぎます。 毎回、c [i * n + j]に追加してメモリにアクセスすることはオプションです。これをすべて変数に追加してから、c [i * n + j]と同等にすることができます。 また、コンパイラーはコードを非常にうまく最適化しますが、パフォーマンスを改善するために分岐点（問題がある場合など）で問題が発生するため、命令パイプラインの深さを増やすようにしてください。 したがって、アルゴリズムの線形部分（計算が条件またはサイクルなしで連続して行われる場合）であるため、3番目のサイクル内で1回ではなく4回の乗算が行われるようにメインサイクルを拡張します。 さて、できるだけ少ないメモリにアクセスして乗算し、ローカル変数とポインターにすべてを保存しようとするので、作業も高速化されます。 ループ内のパリティをチェックするのではなく、一度だけチェックするように、偶数サイズと奇数サイズの行列のコードダビングを行います。 この関数を取得します：
int multMatrixSq（double * a、const int n、double * b、double * c）{
double s1 = 0、s2 = 0、s3 = 0、s4 = 0;
double * f、* g、* h;
int i、j、k;
if（n％2 == 0）{
for（i = 0; i <n-1; i + = 2）
{
f = a + i * n;
for（j = 0; j <n-1; j + = 2）
{
g = b + j;
for（k = 0; k <n; k ++）{
s1 + = f [k] * g [n * k];
s2 + = f [n + k] * g [n * k];
s3 + = f [k] * g [n * k + 1];
s4 + = f [n + k] * g [n * k + 1];
}
h = c + i * n + j;
h [0] = s1; s1 = 0;
h [n] = s2; s2 = 0;
h [1] = s3; s3 = 0;
h [n + 1] = s4; s4 = 0;
}
}
}
その他{
for（i = 0; i <n-1; i + = 2）
{
f = a + i * n;
for（j = 0; j <n-1; j + = 2）
{
g = b + j;
for（k = 0; k <n; k ++）{
s1 + = f [k] * g [n * k];
s2 + = f [n + k] * g [n * k];
s3 + = f [k] * g [n * k + 1];
s4 + = f [n + k] * g [n * k + 1];
}
h = c + i * n + j;
h [0] = s1; s1 = 0;
h [n] = s2; s2 = 0;
h [1] = s3; s3 = 0;
h [n + 1] = s4; s4 = 0;
}
if（j == n-1）{
g = b + j;
for（k = 0; k <n; k ++）{
s1 + = f [k] * g [n * k];
s2 + = f [n + k] * g [n * k];
}
h = c + i * n + j;
h [0] = s1; s1 = 0;
h [n] = s2; s2 = 0;
}
}
if（i == n-1）{
f = a + i * n;
for（j = 0; j <n-1; j + = 2）
{
g = b + j;
for（k = 0; k <n; k ++）{
s1 + = f [k] * g [n * k];
s3 + = f [k] * g [n * k + 1];
}
h = c + i * n + j;
h [0] = s1; s1 = 0;
h [1] = s3; s3 = 0;
}
if（j == n-1）{
g = b + j;
for（k = 0; k <n; k ++）{
s1 + = f [k] * g [n * k];
}
h = c + i * n + j;
h [0] = s1; s1 = 0;
}
}
}
0を返します。
}
さて、最後の仕上げ。 サイクル内でマトリックスを通過するとき、何らかの形でRAMから情報を取得する必要があります（特に列を読み取る場合）。 RAMからキャッシュまで、情報は行に格納されます（つまり、要素を読み取ると（それに続くいくつかの要素がキャッシュに配置されます）、列を通過するたびに、メモリに登る必要があります）。 これを防ぐために、ほとんどの場合、RAMではなくキャッシュを使用する必要がありました。次のようにします。約100 x 100のサブマトリックスにマトリックスを分割します（サイズは関数の引数で指定されます）。 完全に共有されていなくても怖くないです。下から右に長方形のピースを残します。 部分行列をそれぞれA（1,1）、...、A（n、n）と呼びます。 そして、通常の規則に従って、部分行列から行列を乗算します。 答えが通常の乗算と同じであることを数学的に証明することは難しくありませんが、2つの小さな行列を乗算すると、両方がキャッシュに収まり、はるかに高速に乗算されます。 乗算には、「サブマトリックスを取得」および「サブマトリックスを配置して、それが何であったかを追加する」関数が必要です。 行列をゼロにする機能と、矩形行列を乗算する機能。 なぜなら 長方形の乗算はそれほど多くないので、正方形の乗算と同じように最適化しませんでした。 今ではもちろん、長方形の関数に対して1つの良い関数を記述し、正方形の関数のエイリアスを作成する方が正確だったように思えますが、なぜか覚えていませんが、私はしませんでした。 したがって、最終関数を作成します。
int multMatrixBlock（double * a、const int n、double * b、const int m、double * res、double * c、double * d、double * e）{
int n_block = n / m;
int m_small;
int i、j、k、o;
m_small = n-n_block * m;
clearMatrix（res、n * n）;
for（i = 0; i <= n_block; i ++）
{
for（j = 0; j <= n_block; j ++）
{
for（k = 0; k <= n_block; k ++）
{
if（i <n_block && j <n_block && k <n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixSq（c、m、d、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i <n_block && j <n_block && k == n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixRect（c、m、d、m_small、m、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i == n_block && j <n_block && k <n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixRect（c、m_small、d、m、m、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i == n_block && j <n_block && k == n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixRect（c、m_small、d、m_small、m、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i <n_block && j == n_block && k <n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixRect（c、m、d、m、m_small、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i <n_block && j == n_block && k == n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixRect（c、m、d、m_small、m_small、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i == n_block && j == n_block && k <n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixRect（c、m_small、d、m、m_small、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
else if（i == n_block && j == n_block && k == n_block）{
getSubMatrix（a、n、c、m、i、k）;
getSubMatrix（b、n、d、m、k、j）;
multMatrixSq（c、m_small、d、e）;
setSubMatrixAdded（res、n、e、m、i、j）;
}
}
}
}
0を返します。
}
ヒルベルト行列を初期行列として使用します。 すべての要素は1より小さく、乗算するのに非常に便利です（オーバーフローはありません）。
二重式（int i、int j、int n）{
return 1 /（double（i + j + 1））;
}
プログラムの全文はこちらまたはこちら
次のようにコンパイルされたmatrix_test.cppファイル：
g ++ matrix_test.cpp -O3
私のコンピューターでテストした結果は次のとおりです（サブマトリックスのサイズはiMacとToshibaで200であり、Athlonで115よりも最適ですが、サブマトリックスがなくても最適な側面からではありませんでした）。
テストコンピューター：
ご覧のとおり、パフォーマンスの向上は非常に大きなものです（最新のプロセッサで10倍以上の加速）。 これらすべてを1つのプロセスで検討したことを考慮する価値があります。つまり、デュアルコアCore2Duoは役割を果たさなかったということです。 もちろん、将来的には、マルチスレッドまたはmpiを使用してプログラムを並列化し、1.8〜1.9倍の増加を得ることは正しいでしょう。 また、進歩はまだ止まっておらず、5年前に最新のプロセッサがモデルを追い抜いたことも注目に値します。 一般的に、ご覧のとおり、コードはひどいですが、それだけの価値があります。 そして重要なのは行列乗算法そのものではなく、プログラムを最適化するための可能なテクニックです。