算法设计与分析:第三章 分治 3.5Strassen矩阵乘法
2015-08-05 09:50
471 查看
/* Strassen矩阵乘法: A和B的乘积矩阵C中的元素C[i,j]定义为:C[i][j] = k从1到n 累加A[i][k]*B[k][j] 若依此定义来计算A和B的乘积矩阵C,则每计 算C的一个元素C[i][j],需要做n次乘法和n-1次 加法。因此,算出矩阵C的 个元素所需的计算 时间为O(n3) 分治法: 将矩阵A,B和C中的每一矩阵都分成4个大小相等的子矩阵。可以将方程C=AB重写为: [C11 C12] = [A11 A12] [B11 B12] [C21 C22] [A21 A22] [B21 B22] 可以得到 C11 = A11*B11 + A12*B21 C12 = A11*B12 + A12*B22 C21 = A21*B11 + A22*B21 C22 = A21*B12 + A22*B22 为了降低时间复杂度,必须减少乘法的次数: 从8次乘法降为7次,用了7次对于n/2矩阵乘的递归调用和18次n.2矩阵的加减法计算 M1 = A11(B12 - B22) M2 = (A11+A12)*B22 以下代码是转载的,有空的时候再研究 */ #include <iostream> using namespace std; const int N=4; //常量N用来定义矩阵的大小 void main() { void STRASSEN(int n,float A[] ,float B[] ,float C[] ); void input(int n,float p[] ); void output(int n,float C[] ); //函数声明部分 float A ,B ,C ; //定义三个矩阵A,B,C cout<<"现在录入矩阵A :"<<endl<<endl; input(N,A); cout<<endl<<"现在录入矩阵B :"<<endl<<endl; input(N,B); //录入数组 STRASSEN(N,A,B,C); //调用STRASSEN函数计算 output(N,C); //输出计算结果 getchar(); } void input(int n,float p[] ) //矩阵输入函数 { int i,j; for(i=0;i<n;i++) { cout<<"请输入第"<<i+1<<"行"<<endl; for(j=0;j<n;j++) cin>>p[i][j]; } } 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(float A[] ,float B[] ,float C[] ) //按通常的矩阵乘法计算C=AB的子算法(仅做2阶) { int i,j,t; for(i=0;i<2;i++) //计算A*B-->C for(j=0;j<2;j++) { C[i][j]=0; //计算完一个C[i][j],C[i][j]应重新赋值为零 for(t=0;t<2;t++) C[i][j]=C[i][j]+A[i][t]*B[t][j]; } } void MATRIX_ADD(int n,float X[] ,float Y[] ,float Z[] ) //矩阵加法函数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[] ,float Y[] ,float Z[] ) //矩阵减法函数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(int n,float A[] ,float B[] ,float C[] ) //STRASSEN函数(递归) { float A11 ,A12 ,A21 ,A22 ; float B11 ,B12 ,B21 ,B22 ; float C11 ,C12 ,C21 ,C22 ; float M1 ,M2 ,M3 ,M4 ,M5 ,M6 ,M7 ; float AA ,BB ,MM1 ,MM2 ; int i,j;//,x; if (n==2) MATRIX_MULTIPLY(A,B,C);//按通常的矩阵乘法计算C=AB的子算法(仅做2阶) else { 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); STRASSEN(n/2,A11,BB,M1);//M1=A11(B12-B22) MATRIX_ADD(n/2,A11,A12,AA); STRASSEN(n/2,AA,B22,M2);//M2=(A11+A12)B22 MATRIX_ADD(n/2,A21,A22,AA); STRASSEN(n/2,AA,B11,M3);//M3=(A21+A22)B11 MATRIX_SUB(n/2,B21,B11,BB); STRASSEN(n/2,A22,BB,M4);//M4=A22(B21-B11) MATRIX_ADD(n/2,A11,A22,AA); MATRIX_ADD(n/2,B11,B22,BB); STRASSEN(n/2,AA,BB,M5);//M5=(A11+A22)(B11+B22) MATRIX_SUB(n/2,A12,A22,AA); MATRIX_SUB(n/2,B21,B22,BB); STRASSEN(n/2,AA,BB,M6);//M6=(A12-A22)(B21+B22) MATRIX_SUB(n/2,A11,A21,AA); MATRIX_SUB(n/2,B11,B12,BB); STRASSEN(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 } }
相关文章推荐
- Silverlight 2 (beta1)数据操作(2)——使用ASP.NET Web Service进行数据CRUD操作(下)
- [转]在linux下如何判断是否已经安装某个软件?软件安装在哪个目录
- C语言冒泡排序实现
- http get post传参
- HDU--5344--2015 Multi-University Training Contest 5--MZL's xor
- 一个可以下载所有操作系统的网站
- 算法设计与分析:第三章 分治 3.3二进制大整数的乘法
- 算法设计与分析:第三章 分治 3.4多项式乘积的分治方法
- 算法设计与分析:第三章 分治 3.2正整数划分
- 代码整洁之道--注释,格式
- C语言冒泡排序
- c++加密解密算法用java方法替代
- HTTP POST请求报文格式分析与Java实现文件上传
- htons函数详解
- 6种常见的数据加载模式设计
- double my_atof(char *str)
- 机器学习05(SVM理论推导)
- Codeforces 443B Kuriyama Mirai's Stones
- poj 2065 高斯消元取模解方程组
- android应用内切换语言(无需退出应用)