辅助函数
为了便于代码的编写，以及在函数参数中传二维数组的方便，写了以下三个辅助函数。
//数组寻址辅助函数
int& GetArrayVal( int* pMatrix, int nCol, int i, int j )
{
return *( pMatrix + i * nCol + j);
}
//创建矩阵
void CreateMatrix( int** pMatrix, int nRow, int nCol )
{
*pMatrix = new int[nRow*nCol];
memset( *pMatrix, 0, sizeof(int*)*nRow*nCol );
for( int i = 0; i < nRow; ++i )
{
for( int j = 0; j < nCol; ++j )
{
GetArrayVal( *pMatrix, nCol, i, j ) = rand() % 5;
}
}
}
//销毁矩阵
void DeleteMatrix( int** pMatrix )
{
if ( NULL != *pMatrix )
{
delete *pMatrix;
*pMatrix = NULL;
}
}
普通矩阵乘法
根据矩阵乘法的定义，可以迅速写出程序代码，如下
//常规数组乘法
bool MatrixMuiltiplyGeneral( int* pMatrix1, int nRow1, int nCol1,
int* pMatrix2, int nRow2,int nCol2, int** pResult )
{
if ( nCol1 != nRow2 ) return false;
*pResult = new int[nRow1*nCol2];
for( int i = 0; i < nRow1; ++i )
{
for( int j = 0; j < nCol2; ++j )
{
GetArrayVal( *pResult, nCol2, i, j ) = 0;
for( int k = 0; k < nCol1; ++k )
{
GetArrayVal( *pResult, nCol2, i, j ) += GetArrayVal( pMatrix1, nCol1, i, k)
* GetArrayVal( pMatrix2, nCol2, k, j );
}
}
}
return true;
}
递归方法（矩阵方阵，即n*n）----8次乘法
将矩阵A，B，C中每一个矩阵分块成4个大小相等的子矩阵，每个子矩阵都为（n/2）*(n/2)。由此可将方程C = AB重写为
由此可得
///
//矩阵乘法，递归调用（8次乘法）
/*****************************************************************
* 时 间: [2015年6月22日]
* 作 者：shanql
* 函数描述：数组乘法
* 函数参数：pMatrix1--矩阵1
nLeftIndex1，nTopIndex1--矩阵1左上角索引点（相对于源矩阵pMatrixl）
pMatrix2--矩阵2
nLeftIndex2, nTopIndex2--矩阵2左上角索引点（相对于源矩阵pMatrix2）
nCount--方阵n
nTotalCol--使用矩阵列大小（指针取值使用）
pResult--相乘结果保存
* 函数返回：
*****************************************************************/
bool MatrixMuiltiplyGeneralEx_8( int* pMatrix1, int nLeftIndex1, int nTopIndex1,
int* pMatrix2, int nLeftIndex2, int nTopIndex2,
int nCount, int nTotalCol, int** pResult )
{
*pResult = new int[nCount*nCount];
for( int i = 0; i < nCount; ++i )
{
for( int j = 0; j < nCount; ++j )
{
GetArrayVal( *pResult, nCount, i, j ) = 0;
for( int k = 0; k < nCount; ++k )
{
GetArrayVal( *pResult, nCount, i, j ) += GetArrayVal( pMatrix1, nTotalCol, nLeftIndex1+i, nTopIndex1+k)
* GetArrayVal( pMatrix2, nTotalCol, nLeftIndex2+k, nTopIndex2+j );
}
}
}
return true;
}
//矩阵加法(n*n),结果保存在第一个矩阵中
void MatrixAddOrSub( int* pMatrix1, int* pMatrix2, int nCount, bool bAdd = true )
{
for( int i = 0; i < nCount; ++i )
{
for( int j = 0; j < nCount; ++j )
{
if ( bAdd )
{
GetArrayVal( pMatrix1, nCount, i, j ) += GetArrayVal( pMatrix2, nCount, i, j );
}
else
{
GetArrayVal( pMatrix1, nCount, i, j ) -= GetArrayVal( pMatrix2, nCount, i, j );
}
}
}
}
//递归求解n*n方矩阵之积(8次乘法)
void RecursiveMuiltiplyMatrix( int* pMatrix1, int nLeftIndex1, int nTopIndex1,
int* pMatrix2, int nLeftIndex2, int nTopIndex2,
int nCount, int nTotalCol, int** pResult )
{
if ( nCount == 2 )
{
MatrixMuiltiplyGeneralEx_8( pMatrix1, nLeftIndex1, nTopIndex1,
pMatrix2, nLeftIndex2, nTopIndex2, nCount, nTotalCol, pResult );
}
else
{
//拆分成4个大小相等的子矩阵
int* pResult1 = NULL;
int* pResult2 = NULL;
int* pResult3 = NULL;
int* pResult4 = NULL;
int* pResult5 = NULL;
int* pResult6 = NULL;
int* pResult7 = NULL;
int* pResult8 = NULL;
//A[1][1]*B[1][1];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1, nTopIndex1,
pMatrix2, nLeftIndex2, nTopIndex2, nCount / 2, nTotalCol, &pResult1 );
//A[1][2]*B[2][1];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1, nTopIndex1+nCount/2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2, nCount / 2, nTotalCol, &pResult2 );
//A[1][1]*B[1][2];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1, nTopIndex1,
pMatrix2, nLeftIndex2, nTopIndex2+nCount/2, nCount / 2, nTotalCol, &pResult3 );
//A[1][2]*B[2][2];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1, nTopIndex1+nCount/2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2+nCount/2, nCount / 2, nTotalCol, &pResult4 );
//A[2][1]*B[1][1];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1+nCount/2, nTopIndex1,
pMatrix2, nLeftIndex2, nTopIndex2, nCount / 2, nTotalCol, &pResult5 );
//A[2][2]*B[2][1];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1+nCount/2, nTopIndex1+nCount/2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2, nCount / 2, nTotalCol, &pResult6 );
//A[2][1]*B[1][2];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1+nCount/2, nTopIndex1,
pMatrix2, nLeftIndex2, nTopIndex2+nCount/2, nCount / 2, nTotalCol, &pResult7 );
//A[2][2]*B[2][2];
RecursiveMuiltiplyMatrix( pMatrix1, nLeftIndex1+nCount/2, nTopIndex1+nCount/2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2+nCount/2, nCount / 2, nTotalCol, &pResult8 );
//加法运算
MatrixAddOrSub( pResult1, pResult2, nCount/2 );
MatrixAddOrSub( pResult3, pResult4, nCount/2);
MatrixAddOrSub( pResult5, pResult6, nCount/2 );
MatrixAddOrSub( pResult7, pResult8, nCount/2 );
//构造结果
*pResult = new int[nCount*nCount];
for( int i = 0; i < nCount/2; ++i )
{
for( int j = 0; j <nCount/2; ++j )
{
GetArrayVal( *pResult, nCount, i, j ) = GetArrayVal( pResult1, nCount/2, i, j );
GetArrayVal( *pResult, nCount, i, j+nCount/2 ) = GetArrayVal( pResult3, nCount/2, i, j );
GetArrayVal( *pResult, nCount, i+nCount/2, j ) = GetArrayVal( pResult5, nCount/2, i, j );
GetArrayVal( *pResult, nCount, i+nCount/2, j+nCount/2 ) = GetArrayVal( pResult7, nCount/2, i, j );
}
}
DeleteMatrix( &pResult1 );
DeleteMatrix( &pResult2 );
DeleteMatrix( &pResult3 );
DeleteMatrix( &pResult4 );
DeleteMatrix( &pResult5 );
DeleteMatrix( &pResult6 );
DeleteMatrix( &pResult7 );
DeleteMatrix( &pResult8 );
}
}
递归方法（矩阵方阵，即n*n）----7次乘法
7次乘法运算如下：
做了7次乘法后，再做若干次加、减法即可得到结果
/
///
//7次乘法
/*****************************************************************
* 时 间: [2015年6月22日]
* 作 者：shanql
* 函数描述：数组加减法（n*n方阵）
* 函数参数：pMatrix1--矩阵1
nLeftIndex1，nTopIndex1--矩阵1左上角索引点（相对于源矩阵pMatrixl）
nTotalCol1--当前传入矩阵实际使用的列数
pMatrix2--矩阵2
nLeftIndex2, nTopIndex2--矩阵2左上角索引点（相对于源矩阵pMatrix2）
nTotalCol2--当前传入矩阵实际使用的列数
nCount--当前方阵n
pResult--相乘结果保存
bAdd--加减标记
* 函数返回：
*****************************************************************/
void MatrixAddOrSubEx(int* pMatrix1, int nLeftIndex1, int nTopIndex1, int nTotalCol1,
int* pMatrix2, int nLeftIndex2, int nTopIndex2, int nTotalCol2,
int nCount, int** pResult, bool bAdd )
{
*pResult = new int[nCount*nCount];
for( int i = 0; i < nCount; ++i )
{
for( int j = 0; j < nCount; ++j )
{
if ( bAdd )
{
GetArrayVal( *pResult, nCount, i, j ) = GetArrayVal( pMatrix1, nTotalCol1, nLeftIndex1+i, nTopIndex1+j)
+ GetArrayVal( pMatrix2, nTotalCol2, nLeftIndex2+i, nTopIndex2+j );
}
else
{
GetArrayVal( *pResult, nCount, i, j ) = GetArrayVal( pMatrix1, nTotalCol1, nLeftIndex1+i, nTopIndex1+j)
- GetArrayVal( pMatrix2, nTotalCol2, nLeftIndex2+i, nTopIndex2+j );
}
}
}
}
//常规数组乘法
bool MatrixMuiltiplyGeneralEx_7( int* pMatrix1, int nLeftIndex1, int nTopIndex1, int nTotalCol1,
int* pMatrix2, int nLeftIndex2, int nTopIndex2, int nTotalCol2,
int nCount, int** pResult )
{
*pResult = new int[nCount*nCount];
for( int i = 0; i < nCount; ++i )
{
for( int j = 0; j < nCount; ++j )
{
GetArrayVal( *pResult, nCount, i, j ) = 0;
for( int k = 0; k < nCount; ++k )
{
GetArrayVal( *pResult, nCount, i, j ) += GetArrayVal( pMatrix1, nTotalCol1, nLeftIndex1+i, nTopIndex1+k)
* GetArrayVal( pMatrix2, nTotalCol2, nLeftIndex2+k, nTopIndex2+j );
}
}
}
return true;
}
//递归求解n*n方矩阵之积(7次乘法)
void RecursiveMuiltiplyMatrixEx( int* pMatrix1, int nLeftIndex1, int nTopIndex1,int nTotalCol1,
int* pMatrix2, int nLeftIndex2, int nTopIndex2, int nTotalCol2,
int nCount, int** pResult )
{
if ( nCount == 2 )
{
MatrixMuiltiplyGeneralEx_7( pMatrix1, nLeftIndex1, nTopIndex1, nTotalCol1,
pMatrix2, nLeftIndex2, nTopIndex2, nTotalCol2, nCount, pResult );
}
else
{
//拆分成4个大小相等的子矩阵
int* pResultM1 = NULL;
int* pResultM2 = NULL;
int* pResultM3 = NULL;
int* pResultM4 = NULL;
int* pResultM5 = NULL;
int* pResultM6 = NULL;
int* pResultM7 = NULL;
//M1 = A[1][1]*(B[1][2] - B[2][2])
int* pB12_B22 = NULL;
MatrixAddOrSubEx( pMatrix2, nLeftIndex2, nTopIndex2+nCount/2, nTotalCol2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2+nCount/2, nTotalCol2, nCount/2, &pB12_B22, false );
RecursiveMuiltiplyMatrixEx( pMatrix1, nLeftIndex1, nTopIndex1, nTotalCol1,
pB12_B22, 0, 0, nCount/2, nCount/2, &pResultM1 );
//M2 = (A[1][1] + A[1][2]) * B[2][2];
int* pA11_A12 = NULL;
MatrixAddOrSubEx( pMatrix1, nLeftIndex1, nTopIndex1, nTotalCol1,
pMatrix1, nLeftIndex1, nTopIndex1+nCount/2, nTotalCol1, nCount/2, &pA11_A12, true );
RecursiveMuiltiplyMatrixEx( pA11_A12, 0, 0, nCount/2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2+nCount/2, nTotalCol2, nCount/2, &pResultM2 );
//M3 = (A[2][1] + A[2][2]) * B[1][1];
int* pA21_A22 = NULL;
MatrixAddOrSubEx( pMatrix1, nLeftIndex1+nCount/2, nTopIndex1, nTotalCol1,
pMatrix1, nLeftIndex1+nCount/2, nTopIndex1+nCount/2, nTotalCol1, nCount/2, &pA21_A22, true );
RecursiveMuiltiplyMatrixEx( pA21_A22, 0, 0, nCount/2,
pMatrix2, nLeftIndex2, nTopIndex2, nTotalCol2, nCount/2, &pResultM3 );
//M4 = A[2][2] * (B[2][1] - B[1][1])
int* pB21_B11 = NULL;
MatrixAddOrSubEx( pMatrix2, nLeftIndex2+nCount/2, nTopIndex2, nTotalCol2,
pMatrix2, nLeftIndex2, nTopIndex2, nTotalCol2, nCount/2, &pB21_B11, false );
RecursiveMuiltiplyMatrixEx( pMatrix1, nLeftIndex1+nCount/2, nTopIndex1+nCount/2, nTotalCol1,
pB21_B11, 0, 0, nCount/2, nCount/2, &pResultM4 );
//M5 = (A[1][1] + A[2][2])*(B[1][1]+B[2][2])
int* pA11_A22 = NULL;
int* pB11_B22 = NULL;
MatrixAddOrSubEx( pMatrix1, nLeftIndex1, nTopIndex1, nTotalCol1,
pMatrix1, nLeftIndex1+nCount/2, nTopIndex1+nCount/2, nTotalCol1, nCount/2, &pA11_A22, true );
MatrixAddOrSubEx( pMatrix2, nLeftIndex2, nTopIndex2, nTotalCol2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2+nCount/2, nTotalCol2, nCount/2, &pB11_B22, true );
RecursiveMuiltiplyMatrixEx( pA11_A22, 0, 0, nCount/2,
pB11_B22, 0, 0, nCount/2, nCount/2, &pResultM5 );
//M6 = (A[1][2] - A[2][2])*(B[2][1]+B[2][2])
int* pA12_A22 = NULL;
int* pB21_B22 = NULL;
MatrixAddOrSubEx( pMatrix1, nLeftIndex1, nTopIndex1+nCount/2, nTotalCol1,
pMatrix1, nLeftIndex1+nCount/2, nTopIndex1+nCount/2, nTotalCol1, nCount/2, &pA12_A22, false );
MatrixAddOrSubEx( pMatrix2, nLeftIndex2+nCount/2, nTopIndex2, nTotalCol2,
pMatrix2, nLeftIndex2+nCount/2, nTopIndex2+nCount/2, nTotalCol2, nCount/2, &pB21_B22, true );
RecursiveMuiltiplyMatrixEx( pA12_A22, 0, 0, nCount/2,
pB21_B22, 0, 0, nCount/2, nCount/2, &pResultM6 );
//M7 = (A[1][1] - A[2][1])*(B[1][1]+B[1][2])
int* pA11_A21 = NULL;
int* pB11_B12 = NULL;
MatrixAddOrSubEx( pMatrix1, nLeftIndex1, nTopIndex1, nTotalCol1,
pMatrix1, nLeftIndex1+nCount/2, nTopIndex1, nTotalCol1, nCount/2, &pA11_A21, false );
MatrixAddOrSubEx( pMatrix2, nLeftIndex2, nTopIndex2, nTotalCol2,
pMatrix2, nLeftIndex2, nTopIndex2+nCount/2, nTotalCol2, nCount/2, &pB11_B12, true );
RecursiveMuiltiplyMatrixEx( pA11_A21, 0, 0, nCount/2,
pB11_B12, 0, 0, nCount/2, nCount/2, &pResultM7 );
int* pResultC11 = NULL;
int* pResultC12 = NULL;
int* pResultC21 = NULL;
int* pResultC22 = NULL;
int* pResultTemp1 = NULL;
int* pResultTemp2 = NULL;
//C11 = M5 + M4 - M2 + M6
MatrixAddOrSubEx( pResultM5, 0, 0, nCount/2,
pResultM4, 0, 0, nCount/2, nCount/2, &pResultTemp1, true );
MatrixAddOrSubEx( pResultTemp1, 0, 0, nCount/2,
pResultM2, 0, 0, nCount/2, nCount/2, &pResultTemp2, false );
MatrixAddOrSubEx( pResultTemp2, 0, 0, nCount/2,
pResultM6, 0, 0, nCount/2, nCount/2, &pResultC11, true );
DeleteMatrix( &pResultTemp1 );
DeleteMatrix( &pResultTemp2 );
//C12 = M1 + M2
MatrixAddOrSubEx( pResultM1, 0, 0, nCount/2,
pResultM2, 0, 0, nCount/2, nCount/2, &pResultC12, true );
//C21 = M3 + M4
MatrixAddOrSubEx( pResultM3, 0, 0, nCount/2,
pResultM4, 0, 0, nCount/2, nCount/2, &pResultC21, true );
//C22 = M5 + M1 - M3 - M7
MatrixAddOrSubEx( pResultM5, 0, 0, nCount/2,
pResultM1, 0, 0, nCount/2, nCount/2, &pResultTemp1, true );
MatrixAddOrSubEx( pResultTemp1, 0, 0, nCount/2,
pResultM3, 0, 0, nCount/2, nCount/2, &pResultTemp2, false );
MatrixAddOrSubEx( pResultTemp2, 0, 0, nCount/2,
pResultM7, 0, 0, nCount/2, nCount/2, &pResultC22, false );
DeleteMatrix( &pResultTemp1 );
DeleteMatrix( &pResultTemp2 );
//构造结果
*pResult = new int[nCount*nCount];
for( int i = 0; i < nCount/2; ++i )
{
for( int j = 0; j <nCount/2; ++j )
{
GetArrayVal( *pResult, nCount, i, j ) = GetArrayVal( pResultC11, nCount/2, i, j );
GetArrayVal( *pResult, nCount, i, j+nCount/2 ) = GetArrayVal( pResultC12, nCount/2, i, j );
GetArrayVal( *pResult, nCount, i+nCount/2, j ) = GetArrayVal( pResultC21, nCount/2, i, j );
GetArrayVal( *pResult, nCount, i+nCount/2, j+nCount/2 ) = GetArrayVal( pResultC22, nCount/2, i, j );
}
}
//释放内存
DeleteMatrix( &pResultM1 );
DeleteMatrix( &pResultM2 );
DeleteMatrix( &pResultM3 );
DeleteMatrix( &pResultM4 );
DeleteMatrix( &pResultM5 );
DeleteMatrix( &pResultM6 );
DeleteMatrix( &pResultM7 );
// int* pB12_B22 = NULL;
// int* pA11_A12 = NULL;
// int* pA21_A22 = NULL;
// int* pB21_B11 = NULL;
// int* pA11_A22 = NULL;
// int* pB11_B22 = NULL;
// int* pA12_A22 = NULL;
// int* pB21_B22 = NULL;
// int* pA11_A21 = NULL;
// int* pB11_B12 = NULL;
DeleteMatrix( &pB12_B22 );
DeleteMatrix( &pA11_A12 );
DeleteMatrix( &pA21_A22 );
DeleteMatrix( &pB21_B11 );
DeleteMatrix( &pA11_A22 );
DeleteMatrix( &pB11_B22 );
DeleteMatrix( &pA12_A22 );
DeleteMatrix( &pB21_B22 );
DeleteMatrix( &pA11_A21 );
DeleteMatrix( &pB11_B12 );
DeleteMatrix( &pResultTemp1 );
DeleteMatrix( &pResultTemp2 );
DeleteMatrix( &pResultC11 );
DeleteMatrix( &pResultC12 );
DeleteMatrix( &pResultC21 );
DeleteMatrix( &pResultC22 );
}
}
运行结果分析
以上三种求矩阵乘法的方法中，第一种常规解法，易知其算法复杂度为O(n^3)；第二种递归求解（8次相乘），算法复杂度也为O(n^3)；第三种为O(n^2.81)。具体请参见清华大学王晓东的《算法设计与分析》（第2版）的30页----2.5 Strassen矩阵乘法。
理论上第三种方法算法复杂度最小，应该运行起来最快，但是我写出来的代码，并没有验证这个理论。所以，编码还是挺重要的，应该是在递归的过程中不断的new、delete对象，并进行了大量的赋值运算，所以大量的时候消耗在这个上面了。
void PrintMatrix( int* pMatrix, int nRow, int nCol )
{
return;
for( int i = 0; i < nRow; ++i )
{
for( int j = 0; j < nCol; ++j )
{
cout << left << setw( 3 ) << GetArrayVal( pMatrix, nCol, i, j ) << " ";
}
cout << endl;
}
cout << right;
}
int main()
{
clock_t start, finish;
double duration = 0.0;
//srand( (unsigned)time(NULL) );
srand( 0 );
while( true )
{
int* pMatrix1 = NULL;
int* pMatrix2 = NULL;
int* pMatrixResult1 = NULL;
int nRow1, nCol1, nRow2,nCol2;
// nRow1 = 128;
// nCol1 = 128;
// nRow2 = 128;
// nCol2 = 128;
nRow1 = 256;
nCol1 = 256;
nRow2 = 256;
nCol2 = 256;
CreateMatrix( &pMatrix1, nRow1, nCol1 );
CreateMatrix( &pMatrix2, nRow2, nCol2 );
cout << "普通矩阵乘法：";
start = clock();
MatrixMuiltiplyGeneral( pMatrix1, nRow1, nCol1, pMatrix2, nRow2, nCol2, &pMatrixResult1 );
finish = clock();
duration = (double)(finish-start);
cout << "时间消耗:" << duration << "ms" << endl << endl;
DeleteMatrix( &pMatrixResult1 );
pMatrixResult1 = NULL;
cout << "递归矩阵乘法(8次乘法)：";
start = clock();
RecursiveMuiltiplyMatrix( pMatrix1, 0, 0, pMatrix2, 0, 0, nRow1, nRow1, &pMatrixResult1 );
finish = clock();
duration = (double)(finish-start);
cout << "时间消耗:" << duration << "ms" << endl << endl;
DeleteMatrix( &pMatrixResult1 );
pMatrixResult1 = NULL;
cout << "递归矩阵乘法(7次乘法)：";
start = clock();
RecursiveMuiltiplyMatrixEx( pMatrix1, 0, 0, nRow1,
pMatrix2, 0, 0, nRow1, nRow1, &pMatrixResult1 );
finish = clock();
duration = (double)(finish-start);
cout << "时间消耗:" << duration << "ms" << endl << endl;
DeleteMatrix( &pMatrixResult1 );
pMatrixResult1 = NULL;
//MatrixMuiltiplyGeneral( pMatrix1, nRow1, nCol1, pMatrix2, nRow2, nCol2, &pMatrixResult1 );
//RecursiveMuiltiplyMatrix( pMatrix1, 0, 0, pMatrix2, 0, 0, nRow1, nRow1, &pMatrixResult1 );
//RecursiveMuiltiplyMatrixEx( pMatrix1, 0, 0, nRow1, pMatrix2, 0, 0, nRow1, nRow1, &pMatrixResult1 );
PrintMatrix( pMatrix1, nRow1, nCol1 );
cout << "*" << endl;;
PrintMatrix( pMatrix2, nRow2, nCol2);
cout << "=" << endl;
PrintMatrix( pMatrixResult1, nRow1, nCol2 );
//内存释放
DeleteMatrix( &pMatrix1 );
DeleteMatrix( &pMatrix2 );
DeleteMatrix( &pMatrixResult1 );
cout << "-----------------------------------" << endl;
system( "pause" );
}
return 0;
}
作者：山丘儿
