C语言n阶strassen矩阵乘法,算法重拾之路strassen矩阵乘法

出处：http://blog.csdn.net/lttree

第一章：分治与递归

STRASSEN矩阵乘法

算法描述：

矩阵乘法是线性代数中最常见的问题之一，它在数值计算中有广泛的应用。设A 和 B 是两个n × n矩阵，它们的乘积AB同样是一个n×n矩阵。A和B的乘积矩阵C 中的元素cij定义为：

按照这个定义来看，计算A 与 B 矩阵乘法，至少需要 n 次乘法 与 n-1 次加法，所以可以知道，求矩阵C乘法的时间为O(n^3)

算法分析：

strassen矩阵乘法采用类似于 大数乘法 中的分治技术，将计算2个n阶矩阵乘积所需的时间改进到O(n^log7) ≈ O(n^2.81)

这个算法有个前提条件，必须是 两个n×n阶矩阵相乘，而且n必须为2的幂。

这样我们每次都可以把大矩阵分割成四个小矩阵：

由上面图可知：

C11 = A11B11 + A12B21

C12 = A11B12 + A12B22

C21 = A21B11 + A22B21

C22 = A21B12 + A22B22

如果 n=2 ，两个2×2阶矩阵乘法需要 8次乘法 和 4次加法。当子矩阵的阶数大于2时，为求两个子矩阵的乘积，可以继续将矩阵分块，直到子矩阵阶数降为2，这就是分治降阶的递归算法。

按照这个算法，计算2个n阶矩阵的乘积转化为计算8个n/2阶矩阵的乘积和4个 n/2阶 矩阵的加法。 而 2个 n/2阶矩阵的加法显然可以再O(n^2)时间内完成，由此可知上述算法的时间耗费T(n)：

T(n) = O(1) 当n=2

= 8T(n/2)+O(n^2) 当n>2

Ok,解这个递归方程，发现 T(n) = O(n^3)。并没有减少！

这是因为，这个方法并没有减少矩阵的乘法次数，要改进这个算法的时间，必须减少乘法的次数，Strassen就提出了 对于2阶矩阵的乘积方法，只用了7次乘法，用了多次加减，但是效率依然上升了很多。

这 七次乘法为：

M1 = A11( B12 - B22 )

M2 = (A11 + A12)B22

M3 = (A21 +A22)B11

M4 = A22(B21 - B11)

M5 = (A11 + A22)(B11 + B22)

M6 = (A12 - A22)(B21 + B22)

M7 = (A11 - A21)(B11 + B12)

做完这七次乘法，通过一些加减就可以得到，最后矩阵的值：

C11 = M5 + M4 - M2 + M6

C12 = M1 + M2

C21 = M3 + M4

C22 = M5 + M1 - M3 - M7

这样做以后，它的时间T(n)为：

T(n) = O(1) 当n=2

= 7T(n/2)+O(n^2) 当n>2

解得 T(n) = O(n^log7) ≈ O(n^2.81)

程序代码：

const int N = 8; //Define the size of the Matrix

// 输入矩阵

template

void input(int n, T p[][N]) {

for(int i=0; i

for(int j=0; j

cin>>p[i][j];

}

// 输出矩阵

template

void output(int n, T C[][N]) {

for(int i=0; i

for(int j=0; j

cout<

}

cout<

}

}

// 普通的矩阵乘法

template

void Matrix_Mul(T A[][N], T B[][N], T C[][N]) {

for(int i=0; i<2; i++) {

for(int j=0; j<2; j++) {

C[i][j] = 0;

for(int k=0; k<2; k++) {

C[i][j] = C[i][j] + A[i][k]*B[k][j];

}

}

}

}

// 矩阵加法 C=A+B

template

void Matrix_Add(int n, T A[][N], T B[][N], T C[][N]) {

for(int i=0; i

for(int j=0; j

C[i][j] = A[i][j] + B[i][j];

}

}

}

// 矩阵减法 C=A-B

template

void Matrix_Sub(int n, T A[][N], T B[][N], T C[][N]) {

for(int i=0; i

for(int j=0; j

C[i][j] = A[i][j] - B[i][j];

}

}

}

// strassen矩阵乘法

template

void Strassen(int n, T A[][N], T B[][N], T C[][N]) {

// 将第一个矩阵A分割成四个小矩阵

T A11[N][N], A12[N][N], A21[N][N], A22[N][N];

// 将第二个矩阵B分割成四个小矩阵

T B11[N][N], B12[N][N], B21[N][N], B22[N][N];

// M1~M7 先存起来，为后面最终C矩阵做准备

T M1[N][N], M2[N][N], M3[N][N], M4[N][N], M5[N][N], M6[N][N], M7[N][N];

// 四个小矩阵C 最后组合成大矩阵

T C11[N][N], C12[N][N], C21[N][N], C22[N][N];

// 中间变量

T tempA[N][N], tempB[N][N];

if(n == 2) {

Matrix_Mul(A, B, C);

}

else {

//将矩阵A和B分成阶数相同的四个子矩阵，即分治思想。

for(int i=0; i

{

for(int j=0; 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];

}

}

/* ****** 分别计算 M1~M7 ****** */

//Calculate M1 = (A0 + A3) × (B0 + B3)

Matrix_Add(n/2, A11, A22, tempA);

Matrix_Add(n/2, B11, B22, tempB);

Strassen(n/2, AA, BB, M1);

//Calculate M2 = (A2 + A3) × B0

Matrix_Add(n/2, A21, A22, tempA);

Strassen(n/2, tempA, B11, M2);

//Calculate M3 = A0 × (B1 - B3)

Matrix_Sub(n/2, B12, B22, tempB);

Strassen(n/2, A11, tempB, M3);

//Calculate M4 = A3 × (B2 - B0)

Matrix_Sub(n/2, B21, B11, tempB);

Strassen(n/2, A22, tempB, M4);

//Calculate M5 = (A0 + A1) × B3

Matrix_Add(n/2, A11, A12, tempA);

Strassen(n/2, tempA, B22, M5);

//Calculate M6 = (A2 - A0) × (B0 + B1)

Matrix_Sub(n/2, A21, A11, tempA);

Matrix_Add(n/2, B11, B12, tempB);

Strassen(n/2, tempA, tempB, M6);

//Calculate M7 = (A1 - A3) × (B2 + B3)

Matrix_Sub(n/2, A12, A22, tempA);

Matrix_Add(n/2, B21, B22, tempB);

Strassen(n/2, tempA, tempB, M7);

//Calculate C1 = M1 + M4 - M5 + M7

Matrix_Add(n/2, M1, M4, tempA);

Matrix_Sub(n/2, M7, M5, tempB);

Matrix_Add(n/2, tempA, tempB, C11);

//Calculate C2 = M3 + M5

Matrix_Add(n/2, M3, M5, C12);

//Calculate C3 = M2 + M4

Matrix_Add(n/2, M2, M4, C21);

//Calculate C4 = M1 - M2 + M3 + M6

Matrix_Sub(n/2, M1, M2, tempA);

Matrix_Add(n/2, M3, M6, tempB);

Matrix_Add(n/2, tempA, tempB, C22);

//将四个子矩阵组合回大矩阵 C[][N]

for(int i=0; i

{

for(int j=0; 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];

}

}

}

}