矩阵乘法的算法实现 [转载]

一般矩阵乘法算法：

原理：矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的栏数（column）和第二个矩阵的列数（row）相同时才有定义。一般单指矩阵乘积时，指的便是一般矩阵乘积。若A为m×n矩阵，B为n×p矩阵，则他们的乘积AB会是一个m×p矩阵。其乘积矩阵的元素如下面式子得出：

代码如下：

struct mat{int n, m;double data[MAXN][MAXN];
};int mul(mat& c, const mat& a, const mat& b){int i, j, k;if (a.m != b.n)return 0;c.n = a.n;c.m = b.m;for (i = 0; i < c.n; i++)for (j = 0; j < c.m; j++)for (c.data[i][j] = k = 0; k < a.m; k++)c.data[i][j] += a.data[i][k] * b.data[k][j];return 1;
}

时间复杂度：O(n³)

传统分治方法：

C = AB

将A, B, C分成相等大小的方块矩阵：

C₁₁=A₁₁B₁₁+A₁₂B₂₁ (2)

C₁₂=A₁₁B₁₂+A₁₂B₂₂ (3)

C₂₁=A₂₁B₁₁+A₂₂B₂₁ (4)

C₂₂=A₂₁B₁₂+A₂₂B₂₂ (5)

如果n=2，则2个2阶方阵的乘积可以直接用(2)-(5)式计算出来，共需8次乘法和4次加法。当子矩阵的阶大于2时，为求2个子矩阵的积，可以继续将子矩阵分块，直到子矩阵的阶降为2。这样，就产生了一个分治降阶的递归算法。依此算法，计算2个n阶方阵的乘积转化为计算8个n/2阶方阵的乘积和4个n/2阶方阵的加法。2个n/2×n/2矩阵的加法显然可以在c*n2/4时间内完成，这里c是一个常数。因此，上述分治法的计算时间耗费T(n)应该满足：

T(2) = b

T(n) = 8T(n / 2) + cn² (n > 2)

由上式可知：分治法的运用，方阵的乘法的算法效率并没有提高！

分治法的设计思想：将一个难以直接解决的大问题，分割成一些规模较小的相同问题，以便各个击破，分而治之。

任何一个可以用计算机求解的问题所需的计算时间都与其规模有关。问题的规模越小，越容易直接求解，解题所需的计算时间也越少。例如，对于n个元素的排。n=2时，只要作一次比较即可排好序。n=3时序问题，当n=1时，不需任何计算只要作3次比较即可，…。而当n较大时，问题就不那么容易处理了。要想直接解决一个规模较大的问题，有时是相当困难的。

分治法所能解决的问题一般具有以下几个特征：

1.可缩性。问题的规模缩小到一定的程度就可以容易地解决；

2.最有子结构性。问题可以分解为若干个规模较小的相同问题；

3.可合性。利用该问题分解出的子问题的解可以合并为该问题的解；

4.独立性。该问题所分解出的各个子问题是相互独立的，即子问题之间不包含公共的子子问题。

分治思想与递归就像一对孪生兄弟，经常同时应用在算法设计中，并由此产生高效的算法！

Strassen算法

传统方法2个2阶方阵相乘需要8次乘法。按照上述分治法的思想可以看出，要想减少乘法运算次数，关键在于计算2个2阶方阵的乘积时，能否用少于8次的乘法运算。Strassen提出了一种新的算法来计算2个2阶方阵的乘积。他的算法只用了7次乘法运算，但增加了加、减法的运算次数。这7次乘法是:

M₁=A₁₁(B₁₂-B₂₂)

M₂=(A₁₁+A₁₂)B₂₂

M3=(A₂₁+A₂₂)B₁₁

M₄=A₂₂(B₂₁-B₁₁)

M₅=(A₁₁+A₂₂)(B₁₁+B₂₂)

M₆=(A₁₂-A₂₂)(B₂₁+B₂₂)

M₇=(A₁₁-A₂₁)(B₁₁+B₁₂)

做了7次乘法后，再做若干次加、减法：　

C₁₁=M₅+M₄-M₂+M₆

C₁₂=M₁+M₂

C₂₁=M₃+M₄

C₂₂=M₅+M₁-M₃-M₇

算法效率分析:

Strassen矩阵乘积分治算法中，用了7次对于n/2阶矩阵乘积的递归调用和18次n/2阶矩阵的加减运算。由此可知，该算法的所需的计算时间T(n)满足如下的递归方程:

T(2) = b

T(n) = 7T(n / 2) + cn² (n > 2)

按照解递归方程的套用公式法，其解为T(n)=O(n^log7)≈O(n^2.81)。由此可见，Strassen矩阵乘法的计算时间复杂性比普通矩阵乘法有阶的改进。

数据的存储结构：

二维数组（弊端：必须宏定义一常量N）

二重指针（摆脱宏定义的限制，但程序变得繁琐）

原始矩阵:

A,B,C(C=AB)

分块矩阵：

A₁₁，A₁₂，A₂₁，A₂₂(对方阵A分成四块)

B₁₁，B₁₂，B₂₁，B₂₂(对方阵B分成四块)

七块关键阵：

M₁，M₂，M₃，M₄，M₅，M₆，M₇

中间方阵：

AA，BB(计算加法时作一个桥梁)

二级指针定义说明：

int **c;c = (int**)malloc(n * sizeof(int));for(int i = 0; i < n; i++)c[i] = (int *)malloc(n * sizeof(int));

两矩阵相加问题:

//矩阵加法：
void add(int n, int A[][N], int B[][N], int R[][N])
{ int i, j;for(i = 0; i < n; i++)for(j = 0; j < n; j++)R[i][j] = A[i][j] + B[i][j];
}

两矩阵相减问题:

//矩阵减法：
void sub(int n,int A[][N],int B[][N],int R[][N])
{ int i,j;for(i = 0; i < n; i++)for(j = 0; j < n; j++)R[i][j] = A[i][j] - B[i][j];
}

两方阵相乘问题:

//阶为2的方阵相乘：
void multiply(int A[][N],int B[][N],int R[][N])
{   int M1 = A[0][1] * (B[0][1] - B[1][1]);int M2 = (A[0][0] + A[0][1]) * B[1][1];int M3 = (A[1][0] + A[1][1]) * B[0][0];int M4 = A[1][1] * (B[1][0] - B[0][0]);int M5 = (A[0][0] + A[1][1]) * (B[0][0] + B[1][1]);int M6 = (A[0][1] - A[1][1]) * (B[1][0] + B[1][1]);int M7 = (A[0][0] - A[1][0]) * (B[0][0] + B[0][1]);R[0][0] = M5 + M4 - M2 + M6;R[0][1] = M1 + M2;R[1][0] = M3 + M4;R[1][1] = M5 + M1 - M3 - M7;
}

将大的方阵划分的问题

//分块函数
void divide(int A[][N], int n, int X11[][N], int X12[][N], int X21[][N], int X22[][N] )
{int i, j;for(i = 0; i < (n / 2); i++)for(j = 0; j < (n / 2); j++) {X11[i][j] = A[i][j];X12[i][j] = A[i][j + (n / 2)];X21[i][j] = A[i + (n / 2)][j];X22[i][j] = A[i + (n / 2)][j + (n / 2)];}
}

完整代码如下:

#include<stdio.h>
#include<math.h>#define N 4//二阶矩阵相乘
void Matrix_Multiply(int A[][N], int B[][N], int C[][N]) { for(int i = 0; i < 2; i++) {  for(int j = 0; j < 2; j++) {  C[i][j] = 0;        for(int t = 0; t < 2; t++) {  C[i][j] = C[i][j] + A[i][t] * B[t][j];          }    }          }
}  //矩阵加法：
void add(int n, int A[][N], int B[][N], int R[][N])
{ int i, j;for(i = 0; i < n; i++)for(j = 0; j < n; j++)R[i][j] = A[i][j] + B[i][j];
}//矩阵减法：
void sub(int n, int A[][N], int B[][N], int R[][N])
{ int i,j;for(i = 0; i < n; i++)for(j = 0; j < n; j++)R[i][j] = A[i][j] - B[i][j];
}
void strassen(int n, int A[][N], int B[][N], int C[][N])
{int i, j;int A11[N][N], A12[N][N], A21[N][N], A22[N][N];int B11[N][N], B12[N][N], B21[N][N], B22[N][N];int C11[N][N], C12[N][N], C21[N][N], C22[N][N];int AA[N][N], BB[N][N];int M1[N][N], M2[N][N], M3[N][N], M4[N][N], M5[N][N], M6[N][N], M7[N][N];if(n == 2) {Matrix_Multiply(A, B, C);} 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];}}sub(n / 2, B12, B22, BB);strassen(n / 2, A11, BB, M1);add(n / 2, A11, A12, AA);strassen(n / 2, AA, B22, M2);add(n / 2, A21, A22, AA);strassen(n / 2, AA, B11, M3);sub(n / 2, B21, B11, BB);strassen(n / 2, A22, BB, M4);add(n / 2, A11, A22, AA);add(n / 2, B11, B22, BB);strassen(n / 2, AA, BB, M5);sub(n / 2, A12, A22, AA);add(n / 2, B21, B22, BB);strassen(n / 2, AA, BB, M6);sub(n / 2, A11, A21, AA);add(n / 2, B11, B12, BB);strassen(n / 2, AA, BB, M7);//C11 = M5 + M4 - M2 + M6add(n / 2, M5, M4, AA);sub(n / 2, M6, M2, BB);add(n / 2, AA, BB, C11);//C12 = M1 + M2add(n / 2, M1, M2, C12);//C21 = M3 + M4add(n / 2, M3, M4, C21);//C22 = M5 + M1 - M3 - M7sub(n / 2, M5, M3, AA);sub(n / 2, M1, M7, BB);add(n / 2, AA, BB, C22);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];          }          } }
}int main(void)
{int A[N][N], B[N][N], C[N][N];printf("input A: \n");for(int i = 0; i < N; i++)for(int j = 0; j < N; j++)scanf("%d", &A[i][j]);printf("input B: \n");for(int i = 0; i < N; i++)for(int j = 0; j < N; j++)scanf("%d", &B[i][j]);strassen(N, A, B, C);printf("C:\n");for(int i = 0; i < N; i++)for(int j = 0; j < N; j++) {printf("%d ", C[i][j]);if(j > 0 && j % 3 == 0)printf("\n");}return 0;
}

Strassen算法虽然把时间代价从O(n³)降到O(n^2.81)，但是进一步分析，后者的系数比1大很多，当n比较小（如n<45）且矩阵中非零元素较少时，不宜采用此方法。Strassen算法仅当n很大的时候才有价值，因此，可以说这方面的研究更偏重于理论价值。