c++中大矩阵乘法计算的效率问题

假设两个大小相同的方阵需要计算乘法：按照矩阵乘法的规则：

先写一段矩阵初始化代码：

#include <iostream>
#include <cstdlib> 
#include <ctime>
using namespace std;
 
void matrix_print(int **a, int n);
 
int main(int argc, char *argv[])
{
	// 定义数组：
	srand(time(0));
	int matrix_n = 10;
	int numberOfRows = matrix_n;
	int numberOfCols = matrix_n;
	int** mat1 = new int* [numberOfRows];  // a矩阵的行数
	int** mat2 = new int* [numberOfRows]; 
	int** mat3 = new int* [numberOfRows]; 
	for(int i=0; i<numberOfRows; i++)
	{
		mat1[i] = new int[numberOfCols];
		mat2[i] = new int[numberOfCols];
		mat3[i] = new int[numberOfCols];
	} 
	
	// 初始化矩阵  1-10之间的随机数 
	for(int i=0; i<numberOfRows; i++)
	{
		for(int j=0; j<numberOfCols; j++)
		{
			mat1[i][j] = 1 + rand()%(10-1+1);
			mat2[i][j] = 1 + rand()%(10-1+1);
		}
	}
	
	matrix_print(mat1, matrix_n);    // 输出矩阵 
	matrix_print(mat2, matrix_n);    // 输出矩阵 
	//matrix_print(mat3, matrix_n);
	
	return 0;
}
 
 
// 输出矩阵
void matrix_print(int **a, int n)
{
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			cout << a[i][j] << " ";
		}
		cout << endl;
	}
	cout << endl;
} 
 

初步测试：

现在将矩阵乘法的循环的潜逃次序改变一下，改写为一个新的函数，两个函数对比如下：

1.matrix_multiply_ijk版本：

// 计算矩阵乘法 ijk
void matrix_multiply_ijk(int **a, int **b, int **c, int n)  // n表示方阵的阶数 
{
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			int sum = 0;
			for(int k=0; k<n; k++)
			{
				sum += a[i][k]*b[k][j];
			}
			c[i][j] = sum;
		}
	} 
}

2.matrix_multiply_ikj版本：

// 计算矩阵乘法 ikj
void matrix_multiply_ikj(int **a, int **b, int **c, int n)  // n表示方阵的阶数 
{
	for(int i=0; i<n; i++)
	{
		for(int k=0; k<n; k++)
		{
			int sum = 0;
			int j;
			for(j=0; j<n; j++)
			{
				sum += a[i][k]*b[k][j];
			}
			c[i][j] = sum;
		}
	} 
}

性能测量：

3.通过在不同的维度下测试两个函数的运行时间：

c++测试代码运行时间的方法：通过像个的时钟数来测量；

在ijk的方法下：

ijk方法下的测量结果(单位：ms)

	1st	2nd	3rd	4th
n=100	6	3	7	4
n=300	104	90	95	92
n=500	474	492	484	469
n=800	2001	2100	2170	2042
n=1200	10615	9821	9689	9703
n=2000	63677	62797	62698	62551

在ikj的方法下：

ikj方法下的测量结果(单位：ms)

	1st	2nd	3rd	4th
n=100	5	3	3	3
n=300	77	79	81	78
n=500	328	351	341	359
n=800	1435	1413	1423	1453
n=1200	4901	4778	4849	4694
n=2000	21664	21501	21968	21947

结论：

可以看到。在ikj的方案下，矩阵乘法的运算速度较快，而且在矩阵阶数n越大的时候，这种差别越是明显。在计算矩阵乘法的过程中，三层的循环嵌套共有六种排列方式，虽然在每种嵌套方式下，都要执行同样数量的操作，但是花费的时间是不同的。这是因为在不同的嵌套方式下，改变了数据的访问模式，进而改变了缓存未命中的数量。最终影响了运行时间。

关于缓存未命中的简单理解：

简单计算机模型：

L1 一级缓存

L2二级缓存

R 寄存器

ALU算术逻辑单元

在程序开始运行时，数据都位于主存中，需要将参与运算的数据从主存移到寄存器再进行运算。如果需要的数据没有在一级缓存，而是在二级缓存，而需要将数据存二级缓存移动到一级缓存，这称为一级缓存未命中，当需要的数据没有在二级缓存中时，此时为二级缓存未命中，则需要将数据从主存移动到二级缓存，再移动到一级缓存。所以可以通过减少缓存未命中的数量，提高程序的运行效率。计算机会采取相应的策略。

完整代码：

#include <iostream>
#include <cstdlib> 
#include <ctime>
#include <time.h>   // 包含时间测量的函数 
using namespace std;
 
void matrix_multiply_ijk(int **a, int **b, int **c, int n);
void matrix_multiply_ikj(int **a, int **b, int **c, int n);
 
void matrix_print(int **a, int n);
 
int main(int argc, char *argv[])
{
	// 定义数组：
	srand(time(0));
	int matrix_n = 2000;     // 修改矩阵的阶数为不同的值 
	int numberOfRows = matrix_n;
	int numberOfCols = matrix_n;
	int** mat1 = new int* [numberOfRows];  // a矩阵的行数
	int** mat2 = new int* [numberOfRows]; 
	int** mat3 = new int* [numberOfRows]; 
	for(int i=0; i<numberOfRows; i++)
	{
		mat1[i] = new int[numberOfCols];
		mat2[i] = new int[numberOfCols];
		mat3[i] = new int[numberOfCols];
	} 
	
	// 初始化矩阵  1-10之间的随机数 
	for(int i=0; i<numberOfRows; i++)
	{
		for(int j=0; j<numberOfCols; j++)
		{
			mat1[i][j] = 1 + rand()%(10-1+1);
			mat2[i][j] = 1 + rand()%(10-1+1);
		}
	}
	
	//matrix_print(mat1, matrix_n);    // 输出矩阵 
	//matrix_print(mat2, matrix_n);    // 输出矩阵 
	//matrix_print(mat3, matrix_n);
	
	double clocks_PerMills = double(CLOCKS_PER_SEC) / 1000.0;   // 常数，每秒钟包含的时钟数 
	clock_t start_time = clock();                            // 开始的时钟数 
	
	// 选择矩阵乘法方案 
	//matrix_multiply_ijk(mat1, mat2, mat3, matrix_n);        // 矩阵乘法
	matrix_multiply_ikj(mat1, mat2, mat3, matrix_n);
	
	double elapseMills = (clock()-start_time) / clocks_PerMills; 
	cout << "The routine run time: " << elapseMills << "ms" << endl;
	
	cout << "clock_perMils: " << clocks_PerMills << endl; 
	//matrix_print(mat3, matrix_n);
	// matrix_multiply_ikj(mat1, mat2, mat3, matrix_n);
	// matrix_print(mat3, matrix_n);
	// 释放内存 
	for(int i=0; i<numberOfRows; i++)
	{
		delete mat1[i];
		delete mat2[i];
		delete mat3[i];
	}
	delete mat1;
	delete mat2;
	delete mat3;
	return 0;
}
 
// 计算矩阵乘法 ijk
void matrix_multiply_ijk(int **a, int **b, int **c, int n)  // n表示方阵的阶数 
{
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			int sum = 0;
			for(int k=0; k<n; k++)
			{
				sum += a[i][k]*b[k][j];
			}
			c[i][j] = sum;
		}
	} 
}
 
// 计算矩阵乘法 ikj
void matrix_multiply_ikj(int **a, int **b, int **c, int n)  // n表示方阵的阶数 
{
	for(int i=0; i<n; i++)
	{
		for(int k=0; k<n; k++)
		{
			int sum = 0;
			int j;
			for(j=0; j<n; j++)
			{
				sum += a[i][k]*b[k][j];
			}
			c[i][j] = sum;
		}
	} 
}
 
 
// 输出矩阵
void matrix_print(int **a, int n)
{
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			cout << a[i][j] << " ";
		}
		cout << endl;
	}
	cout << endl;
} 

----------------------------------------------------end--------------------------------------------------------