C++ 利用硬件加速矩阵乘法

文章目录

在这里插入图片描述

一、矩阵乘法定义

矩阵 A x × y A_{x \times y} Ax×y 和矩阵 B u × v B_{u \times v} Bu×v 相乘的前提条件是 y = = u y==u y==u ，并且相乘后得到的矩阵为 C x × v C_{x \times v} Cx×v（即 A A A 的行和 B B B 的列构成了矩阵 C C C 的行列）；

二、矩阵类封装

我们用 C++ 封装了一个 n × m n \times m n×m 的矩阵类，用二维数组来存储数据，定义如下：

#define MAXN 1000
#define LL __int64

class Matrix {
private:
	int n, m;
	LL** pkData;
public:
	Matrix() : n(0), m(0) {
		pkData = NULL;
	}
	void Alloc() {
		pkData = new LL *[MAXN];                       // 1)
		for (int i = 0; i < MAXN; ++i) {
			pkData[i] = new LL[MAXN];
		}
	}
	void Dealloc() {
		if (pkData) {
			for (int i = 0; i < MAXN; ++i) {           // 2)
				delete [] pkData[i];
			}
			delete[] pkData;
			pkData = NULL;
		}
	}
};

1） p k D a t a pkData pkData 可以认为是一个二维数组（ p k D a t a [ i ] [ j ] pkData[i][j] pkData[i][j] 就是矩阵第 i 行，第 j 列的数据），之所以这里用了二维指针，是因为当 MAXN 很大时，栈上分配不了这么多空间，容易导致栈溢出，所以通过 new 把空间分配在了堆上；
2）释放空间的时候，首先释放低维空间，再释放高维空间；

三、矩阵乘法实现

1、 i j k ijk ijk 式

最简单的矩阵乘法实现如下：

class Matrix {
	...
public:
	void Multiply_ijk(const Matrix& other, Matrix& ret) {
		// assert(m == other.n);
		ret.Reset(n, other.m);
		int i, j, k;
		for (i = 0; i < n; i++) {
			for (j = 0; j < other.m; j++) {
				for (k = 0; k < m; k++) {
					ret.pkData[i][j] += pkData[i][k] * other.pkData[k][j];
				}
			}
		}
	}
};

这种方法被称为 i j k ijk ijk 式，对矩阵乘法 A × B = C A \times B = C A×B=C，枚举 A A A 的每一行，再枚举 B B B 的每一列，分别对应相乘后放入矩阵 C C C 的对应位置中，如下图所示；

2、 i k j ikj ikj 式

对上述算法进行一些改进，交换两个内层循环的位置，得到如下算法：

class Matrix {
	...
public:
	void Multiply_ikj(const Matrix& other, Matrix& ret) {
		// assert(m == other.n);
		ret.Reset(n, other.m);
		int i, j, k;
		for (i = 0; i < n; i++) {
			for (k = 0; k < m; k++) {
				LL v = pkData[i][k];
				for (j = 0; j < other.m; j++) {
					ret.pkData[i][j] += v * other.pkData[k][j];
				}
			}
		}
	}
};

这种方法被称为 i k j ikj ikj 式，对矩阵乘法 A × B = C A \times B = C A×B=C，行优先枚举 A A A 的每一个格子，再枚举 B B B 的每一行，分别对应相乘后放入矩阵 C C C 的对应位置中，每次相乘得到的 C C C 都是部分积，如下图所示，用绿色的深浅来表示这个值是否已经完整求得；

3、 k i j kij kij 式

对上述算法再进行一些改进，交换两个外层循环的位置，得到如下算法：

class Matrix {
	...
public:
	void Multiply_kij(const Matrix& other, Matrix& ret) {
		// assert(m == other.n);
		ret.Reset(n, other.m);
		int i, j, k;
		for (k = 0; k < m; k++) {
			for (i = 0; i < n; i++) {
				LL v = pkData[i][k];
				for (j = 0; j < other.m; j++) {
					ret.pkData[i][j] += v * other.pkData[k][j];
				}
			}
		}
	}
};

这种方法被称为 k i j kij kij 式，对矩阵乘法 A × B = C A \times B = C A×B=C，列优先枚举 A A A 的每一个格子，再枚举 B B B 的每一行，分别对应相乘后放入矩阵 C C C 的对应位置中，每次相乘得到的 C C C 都是部分积，如下图所示，用绿色的深浅来表示这个值是否已经完整求得；

四、时间测试

矩阵阶数	i j k ijk ijk	i k j ikj ikj	k i j kij kij
200	47 ms	31 ms	16 ms
500	781 ms	438 ms	453 ms
1000	8657 ms	3687 ms	3688 ms
2000	69547 ms	28000 ms	29672 ms

由于矩阵乘法本身的时间复杂度是 O(N3) 的，所以数据量越大，越能看出实际效果；

五、原理分析

原因是因为 CPU 访问内存的速度比 CPU 计算速度慢得多，为了解决速度不匹配的问题，在 CPU 与内存之间加了高速缓存cache。高速缓存 cache 的存在大大提高了 CPU 访问数据的速度。但是当内存访问不连续的时候，就会导致 cache 命中率降低，所以为了加速，就要尽可能使内存访问连续，即不要跳来跳去。
矩阵

六、最后结论

运行速度： i k j ≈ k i j > i j k ikj \approx kij > ijk ikj≈kij>ijk
模板地址：矩阵乘法模板