编程实现矩阵乘法,并考虑当矩阵规模较大时的优化方法:
(1)一般的解法:
//矩阵乘法,3个for循环搞定
void Mul(int** matrixA, int** matrixB, int** matrixC)
{
for(int i = 0; i < n; ++i)
{
for(int j = 0; j < n; ++j)
{
matrixC[i][j] = 0;
for(int k = 0; k < n; ++k)
{
matrixC[i][j] += matrixA[i][k] * matrixB[k][j];
}
}
}
}
(2)strassen算法(分块+分治思想):
算法思想:
strassen算法算法流程:
通用矩阵相乘算法时间复杂度是,而Strassen算法复杂度只是。但随着n的变大,比如当n >> 100时,Strassen算法是比通用矩阵相乘算法变得更有效率。
具体的代码实现如下:
#include <iostream>
#include <ctime>
#include <Windows.h>
using namespace std;
template<typename T>
class Strassen_class{
public:
void ADD(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize); //矩阵相加
void SUB(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize); //矩阵相减
void MUL(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize);//朴素算法实现
void FillMatrix(T** MatrixA, T** MatrixB, int length);//A,B矩阵赋值
void PrintMatrix(T **MatrixA, int MatrixSize);//打印矩阵
void Strassen(int N, T **MatrixA, T **MatrixB, T **MatrixC);//Strassen算法实现
};
//矩阵的相加
template<typename T>
void Strassen_class<T>::ADD(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize)
{
for (int i = 0; i < MatrixSize; i++)
{
for (int j = 0; j < MatrixSize; j++)
{
MatrixResult[i][j] = MatrixA[i][j] + MatrixB[i][j];
}
}
}
// 矩阵的相减
template<typename T>
void Strassen_class<T>::SUB(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize)
{
for (int i = 0; i < MatrixSize; i++)
{
for (int j = 0; j < MatrixSize; j++)
{
MatrixResult[i][j] = MatrixA[i][j] - MatrixB[i][j];
}
}
}
// 矩阵的相乘
template<typename T>
void Strassen_class<T>::MUL(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize)
{
for (int i = 0; i < MatrixSize; i++)
{
for (int j = 0; j < MatrixSize; j++)
{
MatrixResult[i][j] = 0;
for (int k = 0; k < MatrixSize; k++)
{
MatrixResult[i][j] = MatrixResult[i][j] + MatrixA[i][k] * MatrixB[k][j];
}
}
}
}
/*
c++使用二维数组,申请动态内存方法
申请
int **A;
A = new int *[desired_array_row];
for ( int i = 0; i < desired_array_row; i++)
A[i] = new int [desired_column_size];
释放
for ( int i = 0; i < your_array_row; i++)
delete [] A[i];
delete[] A;
*/
template<typename T>
void Strassen_class<T>::Strassen(int N, T **MatrixA, T **MatrixB, T **MatrixC)
{
int HalfSize = N / 2;
int newSize = N / 2;
if (N <= 64) //分治门槛,小于这个值时不再进行递归计算,而是采用常规矩阵计算方法
{
MUL(MatrixA, MatrixB, MatrixC, N);
}
else
{
T** A11;
T** A12;
T** A21;
T** A22;
T** B11;
T** B12;
T** B21;
T** B22;
T** C11;
T** C12;
T** C21;
T** C22;
T** M1;
T** M2;
T** M3;
T** M4;
T** M5;
T** M6;
T** M7;
T** AResult;
T** BResult;
//making a 1 diminsional pointer based array.
A11 = new T *[newSize];
A12 = new T *[newSize];
A21 = new T *[newSize];
A22 = new T *[newSize];
B11 = new T *[newSize];
B12 = new T *[newSize];
B21 = new T *[newSize];
B22 = new T *[newSize];
C11 = new T *[newSize];
C12 = new T *[newSize];
C21 = new T *[newSize];
C22 = new T *[newSize];
M1 = new T *[newSize];
M2 = new T *[newSize];
M3 = new T *[newSize];
M4 = new T *[newSize];
M5 = new T *[newSize];
M6 = new T *[newSize];
M7 = new T *[newSize];
AResult = new T *[newSize];
BResult = new T *[newSize];
int newLength = newSize;
//making that 1 diminsional pointer based array , a 2D pointer based array
for (int i = 0; i < newSize; i++)
{
A11[i] = new T[newLength];
A12[i] = new T[newLength];
A21[i] = new T[newLength];
A22[i] = new T[newLength];
B11[i] = new T[newLength];
B12[i] = new T[newLength];
B21[i] = new T[newLength];
B22[i] = new T[newLength];
C11[i] = new T[newLength];
C12[i] = new T[newLength];
C21[i] = new T[newLength];
C22[i] = new T[newLength];
M1[i] = new T[newLength];
M2[i] = new T[newLength];
M3[i] = new T[newLength];
M4[i] = new T[newLength];
M5[i] = new T[newLength];
M6[i] = new T[newLength];
M7[i] = new T[newLength];
AResult[i] = new T[newLength];
BResult[i] = new T[newLength];
}
//splitting input Matrixes, into 4 submatrices each.
for (int i = 0; i < N / 2; i++)
{
for (int j = 0; j < N / 2; j++)
{
A11[i][j] = MatrixA[i][j];
A12[i][j] = MatrixA[i][j + N / 2];
A21[i][j] = MatrixA[i + N / 2][j];
A22[i][j] = MatrixA[i + N / 2][j + N / 2];
B11[i][j] = MatrixB[i][j];
B12[i][j] = MatrixB[i][j + N / 2];
B21[i][j] = MatrixB[i + N / 2][j];
B22[i][j] = MatrixB[i + N / 2][j + N / 2];
}
}
//here we calculate M1..M7 matrices .
//M1[][]
ADD(A11, A22, AResult, HalfSize);
ADD(B11, B22, BResult, HalfSize); //p5=(a+d)*(e+h)
Strassen(HalfSize, AResult, BResult, M1); //now that we need to multiply this , we use the strassen itself .
//M2[][]
ADD(A21, A22, AResult, HalfSize); //M2=(A21+A22)B11 p3=(c+d)*e
Strassen(HalfSize, AResult, B11, M2); //Mul(AResult,B11,M2);
//M3[][]
SUB(B12, B22, BResult, HalfSize); //M3=A11(B12-B22) p1=a*(f-h)
Strassen(HalfSize, A11, BResult, M3); //Mul(A11,BResult,M3);
//M4[][]
SUB(B21, B11, BResult, HalfSize); //M4=A22(B21-B11) p4=d*(g-e)
Strassen(HalfSize, A22, BResult, M4); //Mul(A22,BResult,M4);
//M5[][]
ADD(A11, A12, AResult, HalfSize); //M5=(A11+A12)B22 p2=(a+b)*h
Strassen(HalfSize, AResult, B22, M5); //Mul(AResult,B22,M5);
//M6[][]
SUB(A21, A11, AResult, HalfSize);
ADD(B11, B12, BResult, HalfSize); //M6=(A21-A11)(B11+B12) p7=(c-a)(e+f)
Strassen(HalfSize, AResult, BResult, M6); //Mul(AResult,BResult,M6);
//M7[][]
SUB(A12, A22, AResult, HalfSize);
ADD(B21, B22, BResult, HalfSize); //M7=(A12-A22)(B21+B22) p6=(b-d)*(g+h)
Strassen(HalfSize, AResult, BResult, M7); //Mul(AResult,BResult,M7);
//C11 = M1 + M4 - M5 + M7;
ADD(M1, M4, AResult, HalfSize);
SUB(M7, M5, BResult, HalfSize);
ADD(AResult, BResult, C11, HalfSize);
//C12 = M3 + M5;
ADD(M3, M5, C12, HalfSize);
//C21 = M2 + M4;
ADD(M2, M4, C21, HalfSize);
//C22 = M1 + M3 - M2 + M6;
ADD(M1, M3, AResult, HalfSize);
SUB(M6, M2, BResult, HalfSize);
ADD(AResult, BResult, C22, HalfSize);
//at this point , we have calculated the c11..c22 matrices, and now we are going to
//put them together and make a unit matrix which would describe our resulting Matrix.
//组合小矩阵到一个大矩阵
for (int i = 0; i < N / 2; i++)
{
for (int j = 0; j < N / 2; j++)
{
MatrixC[i][j] = C11[i][j];
MatrixC[i][j + N / 2] = C12[i][j];
MatrixC[i + N / 2][j] = C21[i][j];
MatrixC[i + N / 2][j + N / 2] = C22[i][j];
}
}
// 释放矩阵内存空间
for (int i = 0; i < newLength; i++)
{
delete[] A11[i]; delete[] A12[i]; delete[] A21[i];
delete[] A22[i];
delete[] B11[i]; delete[] B12[i]; delete[] B21[i];
delete[] B22[i];
delete[] C11[i]; delete[] C12[i]; delete[] C21[i];
delete[] C22[i];
delete[] M1[i]; delete[] M2[i]; delete[] M3[i]; delete[] M4[i];
delete[] M5[i]; delete[] M6[i]; delete[] M7[i];
delete[] AResult[i]; delete[] BResult[i];
}
delete[] A11; delete[] A12; delete[] A21; delete[] A22;
delete[] B11; delete[] B12; delete[] B21; delete[] B22;
delete[] C11; delete[] C12; delete[] C21; delete[] C22;
delete[] M1; delete[] M2; delete[] M3; delete[] M4; delete[] M5;
delete[] M6; delete[] M7;
delete[] AResult;
delete[] BResult;
}//end of else
}
//填充矩阵
template<typename T>
void Strassen_class<T>::FillMatrix(T** MatrixA, T** MatrixB, int length)
{
for (int row = 0; row < length; row++)
{
for (int column = 0; column < length; column++)
{
MatrixB[row][column] = (MatrixA[row][column] = rand() % 5);
}
}
}
//打印矩阵
template<typename T>
void Strassen_class<T>::PrintMatrix(T **MatrixA, int MatrixSize)
{
cout << endl;
for (int row = 0; row < MatrixSize; row++)
{
for (int column = 0; column < MatrixSize; column++)
{
cout << MatrixA[row][column] << "\t";
if ((column + 1) % ((MatrixSize)) == 0)
cout << endl;
}
}
cout << endl;
}
int main()
{
Strassen_class<int> stra;//定义Strassen_class类对象
int MatrixSize = 0;
int** MatrixA; //存放矩阵A
int** MatrixB; //存放矩阵B
int** MatrixC; //存放结果矩阵
clock_t startTime_For_Normal_Multipilication;
clock_t endTime_For_Normal_Multipilication;
clock_t startTime_For_Strassen;
clock_t endTime_For_Strassen;
srand(time(0));
cout << "\n请输入矩阵大小(必须是2的幂指数值(例如:32,64,512,..): ";
cin >> MatrixSize;
int N = MatrixSize;//for readiblity.
//申请内存
MatrixA = new int *[MatrixSize];
MatrixB = new int *[MatrixSize];
MatrixC = new int *[MatrixSize];
for (int i = 0; i < MatrixSize; i++)
{
MatrixA[i] = new int[MatrixSize];
MatrixB[i] = new int[MatrixSize];
MatrixC[i] = new int[MatrixSize];
}
stra.FillMatrix(MatrixA, MatrixB, MatrixSize); //矩阵赋值
//*******************conventional multiplication test
cout << "朴素矩阵算法开始时钟: " << (startTime_For_Normal_Multipilication = clock());
stra.MUL(MatrixA, MatrixB, MatrixC, MatrixSize);//朴素矩阵相乘算法 T(n) = O(n^3)
cout << "\n朴素矩阵算法结束时钟: " << (endTime_For_Normal_Multipilication = clock());
cout << "\n矩阵运算结果... \n";
stra.PrintMatrix(MatrixC, MatrixSize);
//*******************Strassen multiplication test
cout << "\nStrassen算法开始时钟: " << (startTime_For_Strassen = clock());
stra.Strassen(N, MatrixA, MatrixB, MatrixC); //strassen矩阵相乘算法
cout << "\nStrassen算法结束时钟: " << (endTime_For_Strassen = clock());
cout << "\n矩阵运算结果... \n";
stra.PrintMatrix(MatrixC, MatrixSize);
cout << "矩阵大小 " << MatrixSize;
cout << "\n朴素矩阵算法: " << (endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication) << " Clocks.." << (endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication) / CLOCKS_PER_SEC << " Sec";
cout << "\nStrassen算法:" << (endTime_For_Strassen - startTime_For_Strassen) << " Clocks.." << (endTime_For_Strassen - startTime_For_Strassen) / CLOCKS_PER_SEC << " Sec\n";
system("Pause");
return 0;
}
运行结果如下: