strassen算法-现实2次幂,偶数奇数

B12 B22
当矩阵的维数为任一偶数n，总有n=m*2k,将矩阵分为m*m个 的2k阶矩阵。例如n=6,把矩阵分为3*3个21小矩阵。小矩阵用 Strassen相乘，大矩阵用传统算法。
*
=
当矩阵的维数n为奇数时，例如n=5时，分别给矩阵加一行一列 0。
算法代码
指向指针的指针： int **C; C=(int**)malloc(n*sizeof(int)); for(int i=0;i<n;i++) C[i]=(int *)malloc(n*sizeof(int)); 之后，就可以用到C[i][j]取到元素了。
//给矩阵加一行和一列0 void ReinitMatris(){ int i,j; n=n+1; for( i=0;i<n-1;i++){ delete m[i]; } m=new int*[n]; for( i=0;i<n;i++){ m[i]=new int[n]; for( j=0;j<n;j++) m[i][j]=0; } for( i=0;i<n-1;i++) for( j=0;j<n-1;j++) m[i][j]=rand()%10;}
15
//阶为2的n次幂阶矩阵相乘
Matrix Muti2n(Matrix A,Matrix B){ if(A.n==2)//如果阶是2，直接相乘 return MutiFor2(A,B); //将矩阵分成四块 //用Strassen相乘 Matrix M1=Muti2n(Plus(A11,A22),Plus(B11,B22)); Matrix M2=Muti2n(Plus(A21,A22),B11); Matrix M3=Muti2n(A11,Minus(B12,B22)); Matrix M4=Muti2n(A22,Minus(B21,B11)); Matrix M5=Muti2n(Plus(A11,A12),B22); Matrix M6=Muti2n(Minus(A21,A11),Plus(B11,B12)); Matrix M7=Muti2n(Minus(A12,A22),Plus(B21,B22)); Matrix C(A.n); for (int i=0;i<C.n/2;i++) for (int j=0;j<C.n/2;j++){ C.m[i][j]=M1.m[i][j]+M4.m[i][j]-M5.m[i][j]+M7.m[i][j]; C.m[i][j+C.n/2]=M3.m[i][j]+M5.m[i][j]; C.m[i+C.n/2][j]=M2.m[i][j]+M4.m[i][j]; C.m[i+C.n/2][j+C.n/2]=M1.m[i][j]+M3.m[i][j]-M2.m[i][j]+M6.m[i][j]; } return C;}
16
Matrix **cr=new Matrix*[n];//初始化cr，储存ar*br结果 //阶为偶数矩阵想乘 for ( i=0;i<n;i++){ Matrix MutiEven(Matrix A,Matrix B){ cr[i]=new Matrix[n]; int n=A.n,m=1; for ( j=0;j<n;j++) int i,j; cr[i][j]=Matrix(m); //将矩阵阶拆分成m*n，n为奇数 } m为2的幂次 for( i=0;i<n;i++)//用传统方法将ar和br相乘，储存在cr中 while(!(n&1)){ for( j=0;j<n;j++) n>>=1; for(int k=0;k<n;k++) m<<=1; cr[i][j]=Plus(cr[i][j],Muti2n(ar[i][k],br[k][j])); Matrix C(m*n);//声明C } for( i=0;i<n;i++)//将cr合并成一个矩阵C //将A拆分为n^2个m阶矩阵，在ar中储存。 for( j=0;j<n;j++) Matrix **ar=new Matrix*[n]; for(int x=0;x<m;x++) for ( i=0;i<n;i++){ for(int y=0;y<m;y++) ar[i]=new Matrix[n]; for ( j=0;j<n;j++) C.m[i*m+x][j*m+y]=cr[i][j].m[x][y]; ar[i][j]=Split(A,i*m,j*m,m); return C;} } Matrix **br=new Matrix*[n];//同A for (i=0;i<n;i++){ br[i]=new Matrix[n]; for ( j=0;j<n;j++) 17 br[i][j]=Split(B,i*m,j*m,m);}
a13 a23 a33 a43
a14 b11 b12 a24 b21 b22 a34 b31 b32 a44 b41 b42
b13 b23 b33 b43
b14 b24 b34 b44
A11 C A21
A12 B11 A22 B21
6
递归中最了7次（n/2）阶矩阵相乘，18次加法，认为 是一个常数c。 设n阶矩阵相乘需T(n)，则 T(1)=1 T(n)=7T(n/2)+c 解之得：T(n)=O(nlog7)≈O(n2.81)
一般矩阵乘法的时间复杂度为：T(n)=n3。由此可见，Strassen矩阵乘 法的计算时间复杂性比普通矩阵乘法有改进。
7
实现功能：
• 1、实现Strassen算法。
• 2、Strassen要求阶n是2的幂，但这样的情况很少。对任 一偶数n，求n阶矩阵想乘。 • 3、对任意的正整数n，求n阶矩阵相乘。
8
当矩阵的维数2的幂，例如4阶矩阵相乘：C=A*B
a11 a 21 C a31 a41
a12 a22 a32 a42
m1 m2 m3 m4 m5 m6 m7 = = = = = = = (a11 + a22) * (b11 + b22); (a21 + a22) * b11; a11 * (b12 - b22); a22 * (b21 - b11); (a11 + a12) * b22; (a21 - a11) * (b11 + b12); (a12 - a22) * (b21 + b22); 再做若干次加、减法： c11 = m1 + m4 - m5 + m7 c12 = m3 + m5 c21 = m2 + m4 c22 = m1 +m3 -m2 + m6
Leabharlann Baidu
13
//分块函数，分割矩阵A，将第x行第y列的元素作为新矩阵的元素， 分割出一个n阶矩阵。 Matrix Split(Matrix A,int x,int y,int n){ Matrix C(n); for (int i=0;i<n;i++) for (int j=0;j<n;j++) C.m[i][j]=A.m[i+x][j+y]; return C; }
12
//两个矩阵相加，返回和矩阵 //两个矩阵相减，返回差矩阵 Matrix Minus(Matrix A,Matrix B){ Matrix Plus(Matrix A,Matrix B){ Matrix C(A.n); Matrix C(A.n); for (int i=0;i<A.n;i++) for (int i=0;i<A.n;i++) for (int j=0;j<A.n;j++) for (int j=0;j<A.n;j++) C.m[i][j]=A.m[i][j]-B.m[i][j]; C.m[i][j]=A.m[i][j]+B.m[i][j]; return C; return C; } }
Strassen算法
姓名： 王军袖 学号： 31609059
1
目录
1
算法思想
2
3
算法代码
运行结果
2
算法思想
分治算法
其核心思想是将大的问题分解成若干小的子问题进行 处理，从而使复杂的问题变得简单。表现在算法中就 是递归算法。
3
一般的矩阵乘法：
c[i, j ] k 1 a[i, k ] * b[k , j ]
14
//阶为2的方阵相乘： Matrix MutiFor2(Matrix A,Matrix B){//二阶矩阵相乘，返回积 Matrix C(A.n); int M1=(A.m[0][0]+A.m[1][1])*(B.m[0][0]+B.m[1][1]); int M2=(A.m[1][0]+A.m[1][1])*B.m[0][0]; int M3=A.m[0][0]*(B.m[0][1]-B.m[1][1]); int M4=A.m[1][1]*(B.m[1][0]-B.m[0][0]); int M5=(A.m[0][0]+A.m[0][1])*B.m[1][1]; int M6=(A.m[1][0]-A.m[0][0])*(B.m[0][0]+B.m[0][1]); int M7=(A.m[0][1]-A.m[1][1])*(B.m[1][0]+B.m[1][1]); C.m[0][0]=M1+M4-M5+M7; C.m[0][1]=M3+M5; C.m[1][0]=M2+M4; C.m[1][1]=M1+M3-M2+M6; return C;}
k n
C11=A11B11+A12B21 C12=A11B12+A12B22 C21=A21B11+A22B21 C22=A21B12+A22B22
算法分析：
T(n)=n3
传统方法2个2阶方阵相乘需要8次乘法。按照分治的思想可以看 出，要想减少乘法运算次数，关键在于计算2个2阶方阵的乘积时，能 否用少于8次的乘法运算。Strassen提出了一种新的算法来计算2个2阶 方阵的乘积。此算法只用了7次乘法运算，但增加了加、减法的运算 次数。 先计算下面7个量：
5
给出完整的Strassen算法：
procedure STRASSEN(n,A,B,C)(其中n为2的幂) if n=2 then MATRIX-MULTIPLY(A，B，C) else 将矩阵A和B分块; STRASSEN(n/2,A11,B12-B22,M1); STRASSEN(n/2,A11+A12,B22,M2); STRASSEN(n/2,A21+A22,B11,M3); STRASSEN(n/2,A22,B21-B11,M4); STRASSEN(n/2,A11+A22,B11+B22,M5); STRASSEN(n/2,A12-A22,B21+B22,M6); STRASSEN(n/2,A11-21, B11+B12,M7 );
18
当矩阵的维数n是2的幂时，用Multi2n方法，当n为偶数时，用 MutiEven方法，当n为奇数数时，先加一行一列0 ，调用MutiEven方法。 int x=n; while(x%2==0) x=x/2; if(x==1){//2^k阶矩阵 C=Muti2n(A,B); } else if(n%2==0){//偶数阶矩阵 C=MutiEven(A,B); } else{//奇数阶矩阵 A.ReinitMatris(); B.ReinitMatris(); 19 C=MutiEven(A,B);}
运行结果
20
Thank you
23
#include <iostream> #include <time.h> using namespace std; class Matrix{ public: int **m; int n; Matrix(){ m=NULL; n=0; } Matrix(int n){ this->n=n; m=new int*[n]; for(int i=0;i<n;i++){ m[i]=new int[n]; for(int j=0;j<n;j++) m[i][j]=0; } } void InitMatrix(){ for(int i=0;i<n;i++) for(int j=0;j<n;j++) m[i][j]=rand()%10; }