利用SMID指令集对矩阵乘法进行进一步优化

想去卖防晒油的象龟

已于 2023-05-12 09:10:47 修改

阅读量1.2k

点赞数 2

分类专栏：计算机系统原理文章标签：矩阵 c++ 算法

于 2023-05-11 23:31:32 首次发布

本文链接：https://blog.csdn.net/qq_63522621/article/details/130632863

版权

计算机系统原理专栏收录该内容

3 篇文章 0 订阅

订阅专栏

利用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 <windows.h>
#include <iostream>
#include <ctime>
#include <xmmintrin.h>
#include "immintrin.h"
#define N   2048
#define MAX_THREADS 10
using 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 time
def 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矩阵乘法