By SmartPtr(http://www.cppblog.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 阅读(5434)
评论(6) 编辑 收藏 引用