矩阵相乘的快速算法

http://dev.gameres.com/Program/Visual/3D/Mxquick.htm

 

 

赵超 龚敏敏 

 

矩阵相乘的快速算法


算法介绍

矩阵相乘在进行3D变换的时候是经常用到的。在应用中常用矩阵相乘的定义算法对其进行计算。这个算法用到了大量的循环和相乘运算,这使得算法效率不高。而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度,所以对矩阵相乘算法进行一些改进是必要的。

这里要介绍的矩阵算法称为斯特拉森方法,它是由v.斯特拉森在1969年提出的一个方法。

我们先讨论二阶矩阵的计算方法。

对于二阶矩阵

  a11a12    b11b12 
A =a21a22B =b21b22

先计算下面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)来计算出最后结果。

  ma11ma12    mb11mb12 
A4 =ma21ma22B4 =mb21mb22

其中

  a11a12    a13a14    b11b12    b13b14 
ma11 =a21a22ma12 =a23a24mb11 =b21b22mb12 =b23b24
 
  a31a32    a33a34    b31b32    b33b34 
ma21 =a41a42ma22 =a43a44mb21 =b41b42mb22 =b43b44

实现

//  计算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阶矩阵:

 原算法新算法
加法次数4872(48次加法,24次减法)
乘法次数6449
需要额外空间16 * sizeof(float)28 * sizeof(float)

新算法要比原算法多了24次减法运算,少了15次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值