利用AVX-512编程以及循环展开进一步优化
利用AVX-512指令集可同时对512bit数据进行处理，int类型为32bit，因此理论上最多可以将运行速度提升16倍
未优化
DWORD ThreadProc(LPVOID IpParam)//线程函数，用于计算矩阵乘法{ MYDATA* pmd = (MYDATA*)IpParam; int* A = pmd->A, * B = pmd->B, * C = pmd->C; int begin = pmd->begin, end = pmd->end; for (int index = begin; index < end; index++) { int i = index / N, j = index % N; C[i * N + j] = 0; for (int n = 0; n < N; n++) { C[i * N + j] += A[i * N + n] * B[j * N + n]; } } return 0;} 运行结果：
| 矩阵规模 | 1线程 | 2线程 | 4线程 | 8线程 | 16线程 | 32线程 | 64线程 | 100线程 | 1000线程 |
|---|---|---|---|---|---|---|---|---|---|
| 1024*1024 | 2.692s | 1.299s | 0.716s | 0.441s | 0.339s | 0.365s | 0.362s | 0.01s | 0.022s |
| 2048*2048 | 20.407s | 10.593s | 7.006s | 3.879s | 2.725s | 2.787s | 2.836s | 0.012s | 0.017s |
| 4096*4096 | 161.18s | 94.54s | 52.838s | 30.981s | 22.328ss | 22.45s | 22.676s | 0.01s | 0.017s |
矩阵规模
1024
2048
4096
优化后
DWORD ThreadProc(LPVOID IpParam)//线程函数，用于计算矩阵乘法{ MYDATA* pmd = (MYDATA*)IpParam; int* A = pmd->A, * B = pmd->B, * C = pmd->C; int begin = pmd->begin, end = pmd->end; for (int index = begin; index < end; index += 1) { int i = index / N, j = index % N; __m512i res_vec = _mm512_setzero_si512(); for (int n = 0; n < N; n += 16)//循环展开 { __m512i m1_vec = _mm512_loadu_si512((__m512i*) & A[i * N + n]); __m512i m2_vec = _mm512_loadu_si512((__m512i*) & B[j * N + n]); res_vec = _mm512_add_epi32(res_vec, _mm512_mullo_epi32(m1_vec, m2_vec)); } int* p1 = (int*)&res_vec; C[i * N + j] += (p1[0] + p1[1] + p1[2] + p1[3] + p1[4] + p1[5] + p1[6] + p1[7]+p1[8] + p1[9] + p1[10] + p1[11] + p1[12] + p1[13] + p1[14] + p1[15]); } return 0;} 运行结果
| 矩阵规模 | 1线程 | 2线程 | 4线程 | 8线程 | 16线程 | 32线程 | 64线程 | 65线程 | 100线程 | 1000线程 |
|---|---|---|---|---|---|---|---|---|---|---|
| 1024*1024 | 0.327s | 0.18s | 0.108s | 0.088s | 0.063s | 0.074s | 0.085s | 0.013s | 0.028s | 0.409s |
| 2048*2048 | 3.818s | 1.602s | 0.882s | 0.543s | 0.358s | 0.352s | 0.368s | 0.01s | 0.009s | 0.079s |
| 4096*4096 | 23.983s | 13.813s | 7.67s | 4.89s | 3.239s | 3.855s | 3.929s | 0.002s | 0.003s | 0.03s |
矩阵规模
1024
2048
4096
对比发现,采用AVX-512指令集以及循环展开优化后，相比与之前，加速了近10倍，但是当线程数超过64之后，运行速度发生了质变，这时检测矩阵发现并未全部计算完毕，猜测可能为线程数太多，超出库函数所能运行的极限，因此线程被提前销毁，导致未能计算完毕。
发现当线程数为16时运行速度最快，查看cpu规格，cpu线程数为16，因此猜测当程序线程数接近或等于cpu线程数时，运行速度最快。当线程数少时，cpu未全部占用，当线程数多时，操作系统调度可能浪费部分时间。
完整代码
#include#include #include #include #include "immintrin.h"#define N 2048#define MAX_THREADS 10using namespace std;struct MYDATA //用于传递多个参数给线程函数{ int begin, end; int* A, * B, * C;};DWORD ThreadProc(LPVOID IpParam)//线程函数，用于计算矩阵乘法{ MYDATA* pmd = (MYDATA*)IpParam; int* A = pmd->A, * B = pmd->B, * C = pmd->C; int begin = pmd->begin, end = pmd->end; for (int index = begin; index < end; index += 1) { int i = index / N, j = index % N; __m512i res_vec = _mm512_setzero_si512(); for (int n = 0; n < N; n += 16) { __m512i m1_vec = _mm512_loadu_si512((__m512i*) & A[i * N + n]); __m512i m2_vec = _mm512_loadu_si512((__m512i*) & B[j * N + n]); res_vec = _mm512_add_epi32(res_vec, _mm512_mullo_epi32(m1_vec, m2_vec)); } int* p1 = (int*)&res_vec; C[i * N + j] += (p1[0] + p1[1] + p1[2] + p1[3] + p1[4] + p1[5] + p1[6] + p1[7]+p1[8] + p1[9] + p1[10] + p1[11] + p1[12] + p1[13] + p1[14] + p1[15]); } return 0;}int* Mul(int A[], int B[]){ int* M = new int[N * N]; for (int index = 0; index < N * N; index++) { int i = index / N, j = index % N; M[i * N + j] = 0; for (int n = 0; n < N; n++) { M[i * N + j] += A[i * N + n] * B[j * N + n]; } } return M;}int* A=new int[N * N];int* B = new int[N * N];int C[N * N];int* revB = new int[N * N];double rate[MAX_THREADS];int threads[MAX_THREADS] = {1,2,4,8,16,32,64,65,100,1000};int main(){ for (int i = 0; i < N * N; i++) { A[i] = rand() % 2; B[i] = rand() % 2; } for (int index = 0; index < N * N; index++) { int i = index / N, j = index % N; revB[i * N + j] = B[j * N + i]; } for (int i = 0; i < MAX_THREADS; i++) { int m = threads[i]; clock_t start, end; start = clock(); HANDLE hThread[N];//初始化临界区 static MYDATA mydt[N]; int temp = (N * N) / m; for (int i = 0; i < m; i++) { mydt[i].A = A, mydt[i].B = revB, mydt[i].C = C; mydt[i].begin = i * temp, mydt[i].end = i * temp + temp; if (i == m - 1)//最后一个线程 { mydt[i].end = N * N; } hThread[i] = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE)ThreadProc, &mydt[i], 0, NULL);//创建线程 } WaitForMultipleObjects(m, hThread, TRUE, INFINITE); end = clock(); cout << m << "个线程时运行时间为" << (double)(end - start) / CLOCKS_PER_SEC << "s" << endl; }}
python
#单线程实现N*N矩阵乘法import timedef mul(A=[],B=[]): n=len(A) C=[[0]*n for row in range(n)] for m in range(n): for i in range(n): for j in range(n): C[m][i]+=A[m][j]*B[j][i] return C;A= [[1] * 1024 for _ in range(1024)]B= [[1] * 1024 for _ in range(1024)]start=time.time()C=mul(A,B)end=time.time()print("运行时间：",end-start)#多线程实现N*N矩阵乘法 运行时间
| 矩阵规模 | 1024*1024 | 2048*2048 | 4096*4096 |
|---|---|---|---|
| 运行时间 | 96.096s | 763.159s | ∞ infty ∞ |
商业库
| 矩阵规模 | 1024*1024 | 2048*2048 | 4096*4096 |
|---|---|---|---|
| 运行时间 | 0.0339068s | 0.265126s | 2.09942s |
最高加速比（对比python）
| 矩阵规模 | 1024*1024 | 2048*2048 | 4096*4096 |
|---|---|---|---|
| 加速比 | 1525.334 | 2168.065 | ∞ infty ∞ |
通过多线程以及AVX-512指令集合对程序优化后，对比python基准程序最终可实现上千倍优化。
与官方库差距原因：
算法复杂度：矩阵乘法算法的实现可能存在复杂度较高的部分。商业库可能使用了更高效的算法或优化技术来降低计算复杂度，从而提高性能。
并行度：商业库可能使用更高级的并行化技术来充分利用多核处理器。你的优化代码使用了线程来并行执行矩阵乘法，但商业库可能有更细粒度的并行化，利用了更多的硬件资源。
内存访问模式：矩阵乘法的性能受到内存访问的影响。商业库可能通过优化内存访问模式，如使用缓存友好的布局或预取技术，减少了内存访问延迟。在你的优化代码中，对于大型矩阵，内存访问可能不够高效，导致性能下降。
编译器优化：商业库可能使用专门优化的编译器，针对特定硬件和算法进行了深度优化。编译器优化级别、指令选择和代码生成等因素都可以对性能产生影响。
SIMD 和指令集：商业库可能使用了更高级的 SIMD（单指令多数据）指令集来进行向量化操作，从而提高计算性能。你的优化代码中使用的是 AVX-512 SIMD 指令集，但某些商业库可能支持更高级别的 SIMD 指令集，例如 AVX2 或 AVX-512 的高级优化。
性能下降。
编译器优化：商业库可能使用专门优化的编译器，针对特定硬件和算法进行了深度优化。编译器优化级别、指令选择和代码生成等因素都可以对性能产生影响。
SIMD 和指令集：商业库可能使用了更高级的 SIMD（单指令多数据）指令集来进行向量化操作，从而提高计算性能。你的优化代码中使用的是 AVX-512 SIMD 指令集，但某些商业库可能支持更高级别的 SIMD 指令集，例如 AVX2 或 AVX-512 的高级优化。
