【算法导论C++代码】Strassen算法
简单方阵矩乘法
SQUARE-MATRIX-MULTIPLY(A,B) 1 n = A.rows 2 let C be a new n*n natrix 3 for i = 1 to n 4 for j =1 to n 5 cij = 0 6 for k=1 to n 7 cij=cij+aik·bkj 8 return C 一个简单的分治算法 SQUARE-MATRIX-MULTIPLY-RECURSIVE(A,B) 1 n = A.rows 2 let C be a new n*n matrix 3 if n==1 4 c11=a11·b11 5 else partition A,B,and C as in equations (4.9) 6 C11=SQUARE-MATRIX-MULTIPLY-RECURSIVE(A11,B11) +SQUARE-MATRIX-MULTIPLY-RECURSIVE(A12,B21) 7 C12=SQUARE-MATRIX-MULTIPLY-RECURSIVE(A11,B12) +SQUARE-MATRIX-MULTIPLY-RECURSIVE(A12,B22) 8 C21=SQUARE-MATRIX-MULTIPLY-RECURSIVE(A21,B11) +SQUARE-MATRIX-MULTIPLY-RECURSIVE(A22,B21) 9 C22=SQUARE-MATRIX-MULTIPLY-RECURSIVE(A21,B12) +SQUARE-MATRIX-MULTIPLY-RECURSIVE(A22,B22) 10 return C 矩阵乘法的Strassen算法 SQUARE-MATRIX-STRASSEN-RECURSIVE(A,B) 1 n=A.rows 2 let C be a new n*n matrix 3 if n==1 4 c11=a11·b11 5 else partition A,B,and C as in equations(4-9) 6 S1=B12-B22; 7 S2=A11+A12; 8 S3=A21+A22; 9 S4=B21-B11; 10 S5=A11+A22; 11 S6=B11+B22; 12 S7=A12-A22; 13 S8=B21+B22; 14 S9=A11-A21; 15 S10=B11+B12; 16 P1=SQUARE-MATRIX-STRASSEN-RECURSIVE(A11,S1); 17 P2=SQUARE-MATRIX-STRASSEN-RECURSIVE(S2,B22); 18 P3=SQUARE-MATRIX-STRASSEN-RECURSIVE(S3,B11); 19 P4=SQUARE-MATRIX-STRASSEN-RECURSIVE(A22,S4); 20 P5=SQUARE-MATRIX-STRASSEN-RECURSIVE(S5,S6); 21 P6=SQUARE-MATRIX-STRASSEN-RECURSIVE(S7,S8); 22 P7=SQUARE-MATRIX-STRASSEN-RECURSIVE(S9,S10); 23 C11=P5+P4-P2+P6; 24 C12=P1+P2; 25 C21=P3+P4; 26 C22=P5+P1-P3-P7; 26 return C;
/*C++代码。书上给的分解矩阵的做法是用角标计算而不用建立新的对象，不过我并没有想到可以不用新建对象而进行递归的办法，所以这里还是和书上有些不一样的。另外因为演示，所以新建了个类，不过这个类并不稳定，仅作测试时了解功能就好*/
Matrix.h class SquareMatrix { public: SquareMatrix(); SquareMatrix(int **data,int rows); SquareMatrix(int rows); ~SquareMatrix(); int CreateSqMa(int rows); int SetData(int rows,int *data); int **iData; int iRows; friend SquareMatrix operator+(SquareMatrix A,SquareMatrix B); friend SquareMatrix operator-(SquareMatrix A,SquareMatrix B); int SprintSqMa(); }; Matrix.cpp #include <iostream> #include "Matrix.h" SquareMatrix::SquareMatrix() { } SquareMatrix::SquareMatrix(int rows) { this->CreateSqMa(rows); } SquareMatrix::SquareMatrix(int **data,int rows) { iData=data; iRows=rows; } SquareMatrix::~SquareMatrix() { } int SquareMatrix::CreateSqMa(int rows) { iRows = rows; iData =new int *[rows]; for (int i=0;i<iRows;i++) { iData[i]=new int [rows] ; for (int j=0;j<iRows;j++) { iData[i][j]=0; } } return 0; } int SquareMatrix::SetData(int rows,int *data) { int length=rows; for (int i = 0; i < length; i++) { for (int j = 0; j < length; j++) { iData[j][i]=data[i*rows+j]; } } iRows=rows; return 0; } SquareMatrix operator+(SquareMatrix A,SquareMatrix B) { SquareMatrix C(A.iRows); for(int i=0;i<B.iRows;i++) { for(int j=0;j<B.iRows;j++) { C.iData[i][j]=A.iData[i][j]+B.iData[i][j]; } } C.iRows=A.iRows; return C; } SquareMatrix operator-(SquareMatrix A,SquareMatrix B) { SquareMatrix C(A.iRows); for(int i=0;i<B.iRows;i++) { for(int j=0;j<B.iRows;j++) { C.iData[i][j]=A.iData[i][j]-B.iData[i][j]; } } return C; } int SquareMatrix::SprintSqMa() { for(int i=0;i<iRows;i++) { for(int j=0;j<iRows;j++) { std::cout<<iData[i][j]<<' '; if(j==(iRows-1)) { std::cout<<std::endl; } } } return 0; } MAIN.cpp #include <iostream> #include "Matrix.h" using namespace std; SquareMatrix SquareMatrixMultiply(SquareMatrix A,SquareMatrix B); SquareMatrix SquareMatrixMultiplyRecursive(SquareMatrix A,SquareMatrix B); SquareMatrix Strassen(SquareMatrix A,SquareMatrix B); int main() { SquareMatrix A,B,C; B.CreateSqMa(4); C.CreateSqMa(4); A.CreateSqMa(4); int arr[16]={1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; B.SetData(4,arr); A.SetData(4,arr); //C=A+B; //C=SquareMatrixMultiply(A,B); C=SquareMatrixMultiplyRecursive(A,B); //C=Strassen(A,B); cout<<"算法导论4.2矩阵乘法Strassen算法"<<endl; A.SprintSqMa(); cout<<endl; B.SprintSqMa(); cout<<endl; C.SprintSqMa(); cout<<endl; system("pause"); return 0; } SquareMatrix SquareMatrixMultiply(SquareMatrix A,SquareMatrix B) { SquareMatrix C(A.iRows); int n=A.iRows; for (int i=0;i<n;i++) { for(int j=0;j<n;j++) { for(int k=0;k<n;k++) { C.iData[i][j]=C.iData[i][j]+A.iData[i][k]*B.iData[k][j]; } } } return C; } SquareMatrix SquareMatrixMultiplyRecursive(SquareMatrix A,SquareMatrix B) { SquareMatrix C(A.iRows); int n=A.iRows; if(n==1) { C.iData[0][0]=A.iData[0][0]*B.iData[0][0]; } else { int rows_n=n/2; SquareMatrix A11(rows_n),A12(rows_n), A21(rows_n),A22(rows_n), B11(rows_n),B12(rows_n), B21(rows_n),B22(rows_n), C11(rows_n),C12(rows_n), C21(rows_n),C22(rows_n); for (int i=0;i<rows_n;i++) { for(int j=0;j<rows_n;j++) { A11.iData[i][j]=A.iData[i][j]; A12.iData[i][j]=A.iData[i][j+rows_n]; A21.iData[i][j]=A.iData[i+rows_n][j]; A22.iData[i][j]=A.iData[i+rows_n][j+rows_n]; B11.iData[i][j]=B.iData[i][j]; B12.iData[i][j]=B.iData[i][j+rows_n]; B21.iData[i][j]=B.iData[i+rows_n][j]; B22.iData[i][j]=B.iData[i+rows_n][j+rows_n]; } } C11=SquareMatrixMultiplyRecursive(A11,B11) +SquareMatrixMultiplyRecursive(A12,B21); C12=SquareMatrixMultiplyRecursive(A11,B12) +SquareMatrixMultiplyRecursive(A12,B22); C21=SquareMatrixMultiplyRecursive(A21,B11) +SquareMatrixMultiplyRecursive(A22,B21); C22=SquareMatrixMultiplyRecursive(A21,B12) +SquareMatrixMultiplyRecursive(A22,B22); for (int i=0;i<rows_n;i++) { for(int j=0;j<rows_n;j++) { C.iData[i][j]=C11.iData[i][j]; C.iData[i][j+rows_n]=C12.iData[i][j]; C.iData[i+rows_n][j]=C21.iData[i][j]; C.iData[i+rows_n][j+rows_n]=C22.iData[i][j]; } } } return C; } SquareMatrix Strassen(SquareMatrix A,SquareMatrix B) { SquareMatrix C(A.iRows); int n=A.iRows; if(n==1) { C.iData[0][0]=A.iData[0][0]*B.iData[0][0]; } else { int rows_n=n/2; SquareMatrix A11(rows_n),A12(rows_n), A21(rows_n),A22(rows_n), B11(rows_n),B12(rows_n), B21(rows_n),B22(rows_n), C11(rows_n),C12(rows_n), C21(rows_n),C22(rows_n), S1,S2,S3,S4,S5,S6,S7,S8,S9,S10, P1,P2,P3,P4,P5,P6,P7; for (int i=0;i<rows_n;i++) { for(int j=0;j<rows_n;j++) { A11.iData[i][j]=A.iData[i][j]; A12.iData[i][j]=A.iData[i][j+rows_n]; A21.iData[i][j]=A.iData[i+rows_n][j]; A22.iData[i][j]=A.iData[i+rows_n][j+rows_n]; B11.iData[i][j]=B.iData[i][j]; B12.iData[i][j]=B.iData[i][j+rows_n]; B21.iData[i][j]=B.iData[i+rows_n][j]; B22.iData[i][j]=B.iData[i+rows_n][j+rows_n]; } } S1=B12-B22; S2=A11+A12; S3=A21+A22; S4=B21-B11; S5=A11+A22; S6=B11+B22; S7=A12-A22; S8=B21+B22; S9=A11-A21; S10=B11+B12; P1=Strassen(A11,S1); P2=Strassen(S2,B22); P3=Strassen(S3,B11); P4=Strassen(A22,S4); P5=Strassen(S5,S6); P6=Strassen(S7,S8); P7=Strassen(S9,S10); C11=P5+P4-P2+P6; C12=P1+P2; C21=P3+P4; C22=P5+P1-P3-P7; for (int i=0;i<rows_n;i++) { for(int j=0;j<rows_n;j++) { C.iData[i][j]=C11.iData[i][j]; C.iData[i][j+rows_n]=C12.iData[i][j]; C.iData[i+rows_n][j]=C21.iData[i][j]; C.iData[i+rows_n][j+rows_n]=C22.iData[i][j]; } } } return C; }
