矩陣相乘的快速算法
算法介紹
矩陣相乘在進行3D變換的時候是經常用到的。在應用中常用矩陣相乘的定義算法對其進行計算。這個算法用到了大量的循環和相乘運算，這使得算法效率不高。而矩陣相乘的計算效率很大程度上的影響了整個程式的運作速度，是以對矩陣相乘算法進行一些改進是必要的。
這裡要介紹的矩陣算法稱為斯特拉森方法，它是由v.斯特拉森在1969年提出的一個方法。
我們先讨論二階矩陣的計算方法。
對于二階矩陣
a11 a12 b11 b12
A = a21 a22 B = b21 b22
先計算下面7個量(1)
x1 = (a11 + a22) * (b11 + b22);
x2 = (a21 + a22) * b11;
x3 = a11 * (b12 - b22);
x4 = a22 * (b21 - b11);
x5 = (a11 + a12) * b22;
x6 = (a21 - a11) * (b11 + b12);
x7 = (a12 - a22) * (b21 + b22);
再設C = AB。根據矩陣相乘的規則，C的各元素為(2)
c11 = a11 * b11 + a12 * b21
c12 = a11 * b12 + a12 * b22
c21 = a21 * b11 + a22 * b21
c22 = a21 * b12 + a22 * b22
比較(1)(2)，C的各元素可以表示為(3)
c11 = x1 + x4 - x5 + x7
c12 = x3 + x5
c21 = x2 + x4
c22 = x1 + x3 - x2 + x6
根據以上的方法，我們就可以計算4階矩陣了，先将4階矩陣A和B劃分成四塊2階矩陣，分别利用公式計算它們的乘積，再使用(1)(3)來計算出最後結果。
ma11 ma12 mb11 mb12
A4 = ma21 ma22 B4 = mb21 mb22
其中
a11 a12 a13 a14 b11 b12 b13 b14
ma11 = a21 a22 ma12 = a23 a24 mb11 = b21 b22 mb12 = b23 b24
a31 a32 a33 a34 b31 b32 b33 b34
ma21 = a41 a42 ma22 = a43 a44 mb21 = b41 b42 mb22 = b43 b44
實作
// 計算2X2矩陣
void Multiply2X2(float& fOut_11, float& fOut_12, float& fOut_21, float& fOut_22,
float f1_11, float f1_12, float f1_21, float f1_22,
float f2_11, float f2_12, float f2_21, float f2_22)
{
const float x1((f1_11 + f1_22) * (f2_11 + f2_22));
const float x2((f1_21 + f1_22) * f2_11);
const float x3(f1_11 * (f2_12 - f2_22));
const float x4(f1_22 * (f2_21 - f2_11));
const float x5((f1_11 + f1_12) * f2_22);
const float x6((f1_21 - f1_11) * (f2_11 + f2_12));
const float x7((f1_12 - f1_22) * (f2_21 + f2_22));
fOut_11 = x1 + x4 - x5 + x7;
fOut_12 = x3 + x5;
fOut_21 = x2 + x4;
fOut_22 = x1 - x2 + x3 + x6;
}
// 計算4X4矩陣
void Multiply(CLAYMATRIX& mOut, const CLAYMATRIX& m1, const CLAYMATRIX& m2)
{
float fTmp[7][4];
// (ma11 + ma22) * (mb11 + mb22)
Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3],
m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44,
m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44);
// (ma21 + ma22) * mb11
Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3],
m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44,
m2._11, m2._12, m2._21, m2._22);
// ma11 * (mb12 - mb22)
Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3],
m1._11, m1._12, m1._21, m1._22,
m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44);
// ma22 * (mb21 - mb11)
Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3],
m1._33, m1._34, m1._43, m1._44,
m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22);
// (ma11 + ma12) * mb22
Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3],
m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24,
m2._33, m2._34, m2._43, m2._44);
// (ma21 - ma11) * (mb11 + mb12)
Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3],
m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22,
m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24);
// (ma12 - ma22) * (mb21 + mb22)
Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3],
m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44,
m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44);
// 第一塊
mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0];
mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1];
mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2];
mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3];
// 第二塊
mOut._13 = fTmp[2][0] + fTmp[4][0];
mOut._14 = fTmp[2][1] + fTmp[4][1];
mOut._23 = fTmp[2][2] + fTmp[4][2];
mOut._24 = fTmp[2][3] + fTmp[4][3];
// 第三塊
mOut._31 = fTmp[1][0] + fTmp[3][0];
mOut._32 = fTmp[1][1] + fTmp[3][1];
mOut._41 = fTmp[1][2] + fTmp[3][2];
mOut._42 = fTmp[1][3] + fTmp[3][3];
// 第四塊
mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0];
mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1];
mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2];
mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3];
}
比較
在标準的定義算法中我們需要進行n * n * n次乘法運算，新算法中我們需要進行7log2n次乘法，對于最常用的4階矩陣： 原算法 新算法
加法次數 48 72(48次加法，24次減法)
乘法次數 64 49
需要額外空間 16 * sizeof(float) 28 * sizeof(float)
新算法要比原算法多了24次減法運算，少了15次乘法。但因為浮點乘法的運算速度要遠遠慢于加/減法運算，是以新算法的整體速度有所提高。
