By SmartPtr(http://www.shnenglu.com/SmartPtr/)
矩陣相乘在3D變換中是被頻繁用到的一種計算，但在矩陣相乘過程中用到了大量的乘法運算，而cpu中運算單元對于乘法的效率是比較低的，遠低于加法運算，所以，如果能找到一種用加法來替代乘法的方法實現矩陣相乘，將能大大提高我們程序的效率。我們的確有這種方法，這就是網上甚為流行的斯特拉森矩陣乘法，它是由v.斯特拉森在1969年提出的一個方法。
下面對其進行詳細介紹.
一，推導
對于二階矩陣
A = [a11 a12]
[a21 a22]
B = [b11 b12]
[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)來計算出最后結果。
A4 = [ma11 ma12]
[ma21 ma22]
B4 = [mb11 mb12]
[mb21 mb22]
其中
ma11 = [a11 a12]
[a21 a22]
ma12 = [a13 a14]
[a23 a24]
ma21 = [a31 a32]
[a41 a42]
ma22 = [a33 a34]
[a43 a44]
mb11 = [b11 b12]
[b21 b22]
mb12 = [b13 b14]
[b23 b24]
mb21 = [b31 b32]
[b41 b42]
mb22 = [b33 b34]
[b43 b44]
二，實現
typedef float Matrix22[2][2];
typedef float Matrix44[4][4];
inline void Matrix22MulMatrix22(Matrix22 c, const Matrix22& a, const Matrix22& b)
{
float x1 = (a[0][0] + a[1][1]) * (b[0][0] + b[1][1]);
float x2 = (a[1][0] + a[1][1]) * b[0][0];
float x3 = a[0][0] * (b[0][1] - b[1][1]);
float x4 = a[1][1] * (b[1][0] - b[0][0]);
float x5 = (a[0][0] + a[0][1]) * b[1][1];
float x6 = (a[1][0] - a[0][0]) * (b[0][0] + b[0][1]);
float x7 = (a[0][1] - a[1][1]) * (b[1][0] + b[1][1]);
c[0][0] = x1 + x4 -x5 + x7;
c[0][1] = x3 + x5;
c[1][0] = x2 + x4;
c[1][1] = x1 + x3 - x2 + x6;
}
inline void Matrix44MulMatrix44(Matrix44 c, const Matrix44& a, const Matrix44& b)
{
Matrix22 x[7];
// (ma11 + ma22) * (mb11 + mb22)
Matrix22 a0 = {a[0][0]+a[2][2], a[0][1]+a[2][3], a[1][0]+a[3][2], a[1][1]+a[3][3]};
Matrix22 b0 = {b[0][0]+b[2][2], b[0][1]+b[2][3], b[1][0]+b[3][2], b[1][1]+b[3][3]};
Matrix22MulMatrix22(x[0], a0, b0);
// (ma21 + ma22) * mb11
Matrix22 a1 = {a[2][0]+a[2][2], a[2][1]+a[2][3], a[3][0]+a[3][2], a[3][1]+a[3][3]};
Matrix22 b1 = {b[0][0], b[0][1], b[1][0], b[1][1]};
Matrix22MulMatrix22(x[1], a1, b1);
// ma11 * (mb12 - mb22)
Matrix22 a2 = {a[0][0], a[0][1], a[1][0], a[1][1]};
Matrix22 b2 = {b[0][2]-b[2][2], b[0][3]-b[2][3], b[1][2]-b[3][2], b[1][3]-b[3][3]};
Matrix22MulMatrix22(x[2], a2, b2);
// ma22 * (mb21 - mb11)
Matrix22 a3 = {a[2][2], a[2][3], a[3][2], a[3][3]};
Matrix22 b3 = {b[2][0]-b[0][0], b[2][1]-b[0][1], b[3][0]-b[1][0], b[3][1]-b[1][1]};
Matrix22MulMatrix22(x[3], a3, b3);
// (ma11 + ma12) * mb22
Matrix22 a4 = {a[0][0]+a[0][2], a[0][1]+a[0][3], a[1][0]+a[1][2], a[1][1]+a[1][3]};
Matrix22 b4 = {b[2][2], b[2][3], b[3][2], b[3][3]};
Matrix22MulMatrix22(x[4], a4, b4);
// (ma21 - ma11) * (mb11 + mb12)
Matrix22 a5 = {a[2][0]-a[0][0], a[2][1]-a[0][1], a[3][0]-a[1][0], a[3][1]-a[1][1]};
Matrix22 b5 = {b[0][0]+b[0][2], b[0][1]+b[0][3], b[1][0]+b[1][2], b[1][1]+b[1][3]};
Matrix22MulMatrix22(x[5], a5, b5);
// (ma12 - ma22) * (mb21 + mb22)
Matrix22 a6 = {a[0][2]-a[2][2], a[0][3]-a[2][3], a[1][2]-a[3][2], a[1][3]-a[3][3]};
Matrix22 b6 = {b[2][0]+b[2][2], b[2][1]+b[2][3], b[3][0]+b[3][2], b[3][1]+b[3][3]};
Matrix22MulMatrix22(x[6], a6, b6);
// 第一塊
c[0][0] = x[0][0][0] + x[3][0][0] - x[4][0][0] + x[6][0][0];
c[0][1] = x[0][0][1] + x[3][0][1] - x[4][0][1] + x[6][0][1];
c[1][0] = x[0][1][0] + x[3][1][0] - x[4][1][0] + x[6][1][0];
c[1][1] = x[0][1][1] + x[3][1][1] - x[4][1][1] + x[6][1][1];
// 第二塊
c[0][2] = x[2][0][0] + x[4][0][0];
c[0][3] = x[2][0][1] + x[4][0][1];
c[1][2] = x[2][1][0] + x[4][1][0];
c[1][3] = x[2][1][1] + x[4][1][1];
// 第三塊
c[2][0] = x[1][0][0] + x[3][0][0];
c[2][1] = x[1][0][1] + x[3][0][1];
c[3][0] = x[1][1][0] + x[3][1][0];
c[3][1] = x[1][1][1] + x[3][1][1];
// 第四塊
c[2][2] = x[0][0][0] - x[1][0][0] + x[2][0][0] + x[5][0][0];
c[2][3] = x[0][0][1] - x[1][0][1] + x[2][0][1] + x[5][0][1];
c[3][2] = x[0][1][0] - x[1][1][0] + x[2][1][0] + x[5][1][0];
c[3][3] = x[0][1][1] - x[1][1][1] + x[2][1][1] + x[5][1][1];
}
三，分析
在標準的定義算法中我們需要進行n * n * n次乘法運算，新算法中我們需要進行7log2n次乘法，對于最常用的4階矩陣：
原算法 新算法
加法次數 48 72(48次加法，24次減法)
乘法次數 64 49
需要額外空間 16 * sizeof(float) 28 * sizeof(float) （+2 * 4 * 7 * sizeof(float)）
新算法要比原算法多了24次減法運算，少了15次乘法。但因為浮點乘法的運算速度要遠遠慢于加/減法運算，所以新算法的整體速度有所提高。
四，其他
這里列出了按通常公式計算矩陣乘法的函數，以作參考。感謝我的女朋友幫我完成了這兩個函數:)值得一提的是我女朋友是學文科的，從不知道什么是矩陣，當然也沒寫過程序，但在我稍微指點了一下后，等我洗漱完回來，她已經寫好了，經檢查測試通過，把她高興的...
inline void Matrix22MulMatrix22_(Matrix22 c, const Matrix22& a, const Matrix22& b)
{
c[0][0] = a[0][0] * b[0][0] + a[0][1]*b[1][0];
c[0][1] = a[0][0] * b[0][1] + a[0][1]*b[1][1];
c[1][0] = a[1][0] * b[0][0] + a[1][1]*b[1][0];
c[1][1] = a[1][0] * b[0][1] + a[1][1]*b[1][1];
}
inline void Matrix44MulMatrix44_(Matrix44 c, const Matrix44& a, const Matrix44& b)
{
c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0]+a[0][3]*b[3][0];
c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1]+a[0][3]*b[3][1];
c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2]+a[0][3]*b[3][2];
c[0][3] = a[0][0]*b[0][3]+a[0][1]*b[1][3]+a[0][2]*b[2][3]+a[0][3]*b[3][3];
c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0]+a[1][3]*b[3][0];
c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1]+a[1][3]*b[3][1];
c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2]+a[1][3]*b[3][2];
c[1][3] = a[1][0]*b[0][3]+a[1][1]*b[1][3]+a[1][2]*b[2][3]+a[1][3]*b[3][3];
c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0]+a[2][3]*b[3][0];
c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1]+a[2][3]*b[3][1];
c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2]+a[2][3]*b[3][2];
c[2][3] = a[2][0]*b[0][3]+a[2][1]*b[1][3]+a[2][2]*b[2][3]+a[2][3]*b[3][3];
c[3][0] = a[3][0]*b[0][0]+a[3][1]*b[1][0]+a[3][2]*b[2][0]+a[3][3]*b[3][0];
c[3][1] = a[3][0]*b[0][1]+a[3][1]*b[1][1]+a[3][2]*b[2][1]+a[3][3]*b[3][1];
c[3][2] = a[3][0]*b[0][2]+a[3][1]*b[1][2]+a[3][2]*b[2][2]+a[3][3]*b[3][2];
c[3][3] = a[3][0]*b[0][3]+a[3][1]*b[1][3]+a[3][2]*b[2][3]+a[3][3]*b[3][3];
}
當然， 這個用for循環寫出來要簡潔些，但是，這樣更原汁原味:)
posted on 2007-08-26 20:43
SmartPtr 閱讀(5475)
評論(6) 編輯 收藏 引用