【算法导论】矩阵乘法
2015-02-11 11:38
176 查看
离过年都不到十天了，还要等到这周五才能回家，想想也一年没回家了。从寒假开始到现在，已经有二十来天，这期间把2014年总结中的寒假计划也大多数完成了：The Element Of Style的阅读，三门数学课《随机过程》、《工程优化》、《数值分析》的算法实现。回家过年期间肯定不会写博客了，今天一看，这个月只写了三篇，于是乎今天必须再写一篇来完成这个月的基本工作量。言归正传，这篇文章写写选修课《算法设计》作业题中的矩阵乘法的三种方法。
理论公式
C=AB
cij=∑k=1naikbkj
算法实现
分块相乘法
理论公式
A=[A11A21A12A22],B=[B11B21B12B22]
C=AB=[C11C21C12C22]
C11=A11B11+A12B21C12=A11B12+A12B22C21=A21B11+A22B21C22=A21B12+A22B22
算法实现
Strassen法
理论公式
M1=A11(B12−B22)M2=B22(A11+A12)M3=B11(A21+A22)M4=A22(B21−B11)M5=(A11+A22)(B11+B22)M6=(A12−A22)(B21+B22)M7=(A12−A21)(B11+B12)C11=M4+M5+M6−M2C12=M1+M2C21=M3+M4C22=M1+M5−M3−M7
算法实现
完整程序实现
运行结果如下图：
原文：/article/1509033.html
作者：nineheadedbird
矩阵乘法
传统方法理论公式
C=AB
cij=∑k=1naikbkj
算法实现
void TraditionalMethod(float A[][N],float B[][N],float C[]
)//传统方法，三重循环
{
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
C[i][j]=0;//之所以每次调用都清零，是因为前面是循环调用，如果只调用一次就不需要
for(int i=0;i<N;i++)
{
for(int j=0;j<N;j++)
{
for(int k=0;k<N;k++)
{
C[i][j]=C[i][j]+A[i][k]*B[k][j];
}
}
}
}分块相乘法
理论公式
A=[A11A21A12A22],B=[B11B21B12B22]
C=AB=[C11C21C12C22]
C11=A11B11+A12B21C12=A11B12+A12B22C21=A21B11+A22B21C22=A21B12+A22B22
算法实现
void BlockMatrix()//分块矩阵计算
{
for(int i=0;i<N/2;i++)
for(int j=0;j<N/2;j++)
{
A11[i][j]=A[i][j];
A12[i][j]=A[i][j+N/2];
A21[i][j]=A[i+N/2][j];
A22[i][j]=A[i+N/2][j+N/2];
B11[i][j]=B[i][j];
B12[i][j]=B[i][j+N/2];
B21[i][j]=B[i+N/2][j];
B22[i][j]=B[i+N/2][j+N/2];
C11[i][j]=0;
C12[i][j]=0;
C21[i][j]=0;
C22[i][j]=0;
} //将矩阵A和B式分为四块
MATRIX_Multiply(N/2,A11,B11, AA);
MATRIX_Multiply(N/2,A12,B21, BB);
MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z
MATRIX_Multiply(N/2,A11,B12, AA);
MATRIX_Multiply(N/2,A12,B22, BB);
MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z
MATRIX_Multiply(N/2,A21,B11, AA);
MATRIX_Multiply(N/2,A22,B21, BB);
MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z
MATRIX_Multiply(N/2,A21,B12, AA);
MATRIX_Multiply(N/2,A22,B22, BB);
MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z
for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
for(int j=0;j<N/2;j++)
{
C[i][j]=C11[i][j];
C[i][j+N/2]=C12[i][j];
C[i+N/2][j]=C21[i][j];
C[i+N/2][j+N/2]=C22[i][j];
} //计算结果送回C[N][N]
}Strassen法
理论公式
M1=A11(B12−B22)M2=B22(A11+A12)M3=B11(A21+A22)M4=A22(B21−B11)M5=(A11+A22)(B11+B22)M6=(A12−A22)(B21+B22)M7=(A12−A21)(B11+B12)C11=M4+M5+M6−M2C12=M1+M2C21=M3+M4C22=M1+M5−M3−M7
算法实现
void STRASSEN() //STRASSEN函数
{
int i,j;//,x;
for(i=0;i<N/2;i++)
for(j=0;j<N/2;j++)
{
A11[i][j]=A[i][j];
A12[i][j]=A[i][j+N/2];
A21[i][j]=A[i+N/2][j];
A22[i][j]=A[i+N/2][j+N/2];
B11[i][j]=B[i][j];
B12[i][j]=B[i][j+N/2];
B21[i][j]=B[i+N/2][j];
B22[i][j]=B[i+N/2][j+N/2];
} //将矩阵A和B式分为四块
MATRIX_SUB(N/2,B12,B22,BB);
MATRIX_Multiply(N/2,A11,BB,M1);
MATRIX_ADD(N/2,A11,A12,AA);
MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22
MATRIX_ADD(N/2,A21,A22,AA);
MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11
MATRIX_SUB(N/2,B21,B11,BB);
MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)
MATRIX_ADD(N/2,A11,A22,AA);
MATRIX_ADD(N/2,B11,B22,BB);
MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)
MATRIX_SUB(N/2,A12,A22,AA);
MATRIX_ADD(N/2,B21,B22,BB);
MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)
MATRIX_SUB(N/2,A11,A21,AA);
MATRIX_ADD(N/2,B11,B12,BB);
MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
//计算M1,M2,M3,M4,M5,M6,M7（递归部分）
MATRIX_ADD(N/2,M5,M4,MM1);
MATRIX_SUB(N/2,M2,M6,MM2);
MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6
MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2
MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4
MATRIX_ADD(N/2,M5,M1,MM1);
MATRIX_ADD(N/2,M3,M7,MM2);
MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7
for(i=0;i<N/2;i++)
for(j=0;j<N/2;j++)
{
C[i][j]=C11[i][j];
C[i][j+N/2]=C12[i][j];
C[i+N/2][j]=C21[i][j];
C[i+N/2][j+N/2]=C22[i][j];
} //计算结果送回C[N][N]
}完整程序实现
#include <iostream>
#include<ctime>
using namespace std;
const int N=32; //常量N用来定义方阵的大小
void output(int n,float C[]
); //函数声明部分
void TraditionalMethod(float A[][N],float B[][N],float C[]
);//传统的矩阵相乘
void BlockMatrix();
void STRASSEN();
void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2]);
float A[N][N];
float B[N][N];
float C[N][N]; //定义三个矩阵A,B,C
float A11[N/2][N/2],A12[N/2][N/2],A21[N/2][N/2],A22[N/2][N/2];
float B11[N/2][N/2],B12[N/2][N/2],B21[N/2][N/2],B22[N/2][N/2];
float C11[N/2][N/2],C12[N/2][N/2],C21[N/2][N/2],C22[N/2][N/2];
float M1[N/2][N/2],M2[N/2][N/2],M3[N/2][N/2],M4[N/2][N/2],M5[N/2][N/2],M6[N/2][N/2],M7[N/2][N/2];
float AA[N/2][N/2],BB[N/2][N/2],MM1[N/2][N/2],MM2[N/2][N/2];
void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]); //矩阵加法函数X+Y—>Z
void main()
{
//初始化，使相乘的两个矩阵都为全1矩阵
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
{
A[i][j]=1;
B[i][j]=1;
}
//将结果矩阵C初始化为全0矩阵
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
C[i][j]=0;
clock_t start, finish; //用于计时
double duration;
int loop=0;
cout<<"当矩阵为"<<N<<"×"<<N<<"，循环次数为10000时："<<endl<<endl;
//--------使用传统方法--------------//
start = clock();
while(loop<10000)//循环10000次,这里可以更改
{
loop++;
TraditionalMethod(A,B,C); //传统方法计算
}
finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
cout<<"使用传统方法："<<endl;
cout<<"所需时间为："<<duration<<endl<<endl;
// output(N,C); //输出计算结果
//---------使用分块矩阵乘法------------//
start = clock();
loop=0;
while(loop<10000)//循环10000次,这里可以更改
{
loop++;
BlockMatrix(); //分块矩阵计算
}
finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
cout<<"使用分块相乘方法："<<endl;
cout<<"所需时间为："<<duration<<endl<<endl;
// output(N,C); //输出计算结果
//-------使用strassen方法-------------------//
start = clock();
loop=0;
while(loop<10000)//当时间非常小时，需要加大循环次数,这里可以更改
{
loop++;
STRASSEN(); //调用STRASSEN函数计算
}
finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
cout<<"使用strassen方法："<<endl;
cout<<"所需时间为："<<duration<<endl;
// output(N,C); //输出计算结果
}
void TraditionalMethod(float A[][N],float B[][N],float C[] )//传统方法，三重循环 { for(int i=0;i<N;i++) for(int j=0;j<N;j++) C[i][j]=0;//之所以每次调用都清零，是因为前面是循环调用，如果只调用一次就不需要 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { for(int k=0;k<N;k++) { C[i][j]=C[i][j]+A[i][k]*B[k][j]; } } } }
void BlockMatrix()//分块矩阵计算
{
for(int i=0;i<N/2;i++)
for(int j=0;j<N/2;j++)
{
A11[i][j]=A[i][j];
A12[i][j]=A[i][j+N/2];
A21[i][j]=A[i+N/2][j];
A22[i][j]=A[i+N/2][j+N/2];
B11[i][j]=B[i][j];
B12[i][j]=B[i][j+N/2];
B21[i][j]=B[i+N/2][j];
B22[i][j]=B[i+N/2][j+N/2];
C11[i][j]=0;
C12[i][j]=0;
C21[i][j]=0;
C22[i][j]=0;
} //将矩阵A和B式分为四块
MATRIX_Multiply(N/2,A11,B11, AA);
MATRIX_Multiply(N/2,A12,B21, BB);
MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z
MATRIX_Multiply(N/2,A11,B12, AA);
MATRIX_Multiply(N/2,A12,B22, BB);
MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z
MATRIX_Multiply(N/2,A21,B11, AA);
MATRIX_Multiply(N/2,A22,B21, BB);
MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z
MATRIX_Multiply(N/2,A21,B12, AA);
MATRIX_Multiply(N/2,A22,B22, BB);
MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z
for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
for(int j=0;j<N/2;j++)
{
C[i][j]=C11[i][j];
C[i][j+N/2]=C12[i][j];
C[i+N/2][j]=C21[i][j];
C[i+N/2][j+N/2]=C22[i][j];
} //计算结果送回C
}
void output(int n,float C[]
) //矩阵输出函数
{
int i,j;
cout<<"输出矩阵:"<<endl;
for(i=0;i<n;i++)
{
cout<<endl;
for(j=0;j<n;j++)
cout<<C[i][j]<<" ";
}
cout<<endl<<endl;
}
void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2])//矩阵加法函数X*Y—>C
{
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
C[i][j]=0;
for(int k=0;k<n;k++)
{
C[i][j]=C[i][j]+A[i][k]*B[k][j];
}
}
}
}
void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵加法函数X+Y—>Z
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
Z[i][j]=X[i][j]+Y[i][j];
}
void MATRIX_SUB(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵减法函数X-Y—>Z
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
Z[i][j]=X[i][j]-Y[i][j];
}
void STRASSEN() //STRASSEN函数
{
int i,j;//,x;
for(i=0;i<N/2;i++)
for(j=0;j<N/2;j++)
{
A11[i][j]=A[i][j];
A12[i][j]=A[i][j+N/2];
A21[i][j]=A[i+N/2][j];
A22[i][j]=A[i+N/2][j+N/2];
B11[i][j]=B[i][j];
B12[i][j]=B[i][j+N/2];
B21[i][j]=B[i+N/2][j];
B22[i][j]=B[i+N/2][j+N/2];
} //将矩阵A和B式分为四块
MATRIX_SUB(N/2,B12,B22,BB);
MATRIX_Multiply(N/2,A11,BB,M1);
MATRIX_ADD(N/2,A11,A12,AA);
MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22
MATRIX_ADD(N/2,A21,A22,AA);
MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11
MATRIX_SUB(N/2,B21,B11,BB);
MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)
MATRIX_ADD(N/2,A11,A22,AA);
MATRIX_ADD(N/2,B11,B22,BB);
MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)
MATRIX_SUB(N/2,A12,A22,AA);
MATRIX_ADD(N/2,B21,B22,BB);
MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)
MATRIX_SUB(N/2,A11,A21,AA);
MATRIX_ADD(N/2,B11,B12,BB);
MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
//计算M1,M2,M3,M4,M5,M6,M7（递归部分）
MATRIX_ADD(N/2,M5,M4,MM1);
MATRIX_SUB(N/2,M2,M6,MM2);
MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6
MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2
MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4
MATRIX_ADD(N/2,M5,M1,MM1);
MATRIX_ADD(N/2,M3,M7,MM2);
MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7
for(i=0;i<N/2;i++)
for(j=0;j<N/2;j++)
{
C[i][j]=C11[i][j];
C[i][j+N/2]=C12[i][j];
C[i+N/2][j]=C21[i][j];
C[i+N/2][j+N/2]=C22[i][j];
} //计算结果送回C
}
运行结果如下图：
原文：/article/1509033.html
作者：nineheadedbird
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 算法导论之矩阵链乘法详解
- 算法导论15章 动态规划之矩阵链乘法问题
- 【算法导论】矩阵乘法
- 算法提高 矩阵乘法 区间dp
- 算法提高 矩阵乘法
- 蓝桥杯 ALGO-86 算法训练 矩阵乘法
- 算法训练 矩阵乘法
- 算法提高 矩阵乘法 （区间dp）
- 纯C语言矩阵乘法的Strassen算法，包含非2次幂的情况
- 蓝桥杯练习： 算法训练 矩阵乘法
- 51NOD 1836 战忽局的手段（矩阵乘法加速+__float128精度问题+概率期望）——算法马拉松20（告别美国大选及卡斯特罗）
- 转置矩阵的分块并行乘法（C语言实现），计算矩阵C[rawn][rawn]=A[rawm][rawn]'*B[rawm][rawn]，子块大小为S*T，其算法实现原理参加本代码的附件。
- 算法训练 矩阵乘法
- 每日算法--矩阵乘法优化递推
- 算法提高 矩阵乘法 区间DP
- 算法训练 矩阵乘法
- 算法提高 矩阵乘法
- 【数据结构与算法】数组应用1：矩阵乘法（Java实现）
- 程序算法艺术与实践:递归策略之矩阵乘法问题
- （算法学习）矩阵乘法及其应用