高性能计算-gemm-openmp效率测试(10)

1. 目标

设计一个程序，使用OpenMP并行化实现矩阵乘法。给定两个矩阵 A 和 B，矩阵大小均为1024*1024，你的任务是计算它们的乘积 C。
要求：
(1)、使用循环结构体的知识点，包括for循环体并行化、变量规约属性与子句reduction、循环调度策略与子句schedule以及嵌套循环与子句collapse。
(2)、实现OpenMP并行化以加速矩阵乘法的计算。
(3)、考虑内存一致性，确保数据在并行计算过程中的正确性。
(4)、可选：实现线程亲核性，将线程绑定到特定的CPU核心上执行。

2. 测试简介

使用对比串行和并行，及 collapse 和 reduction 的效率。

3. 测试代码

#include <stdio.h>
#include <omp.h>
#include <sys/time.h>
#include <stdlib.h>
 
void gemm_1(int* A,int* B,long* C,int N,int NT);
void gemm_2(int* A,int* B,long* C,int N,int NT);
void gemm_3(int* A,int* B,long* C,int N,int NT);
void gemm_4(int* A,int* B,long* C,int N,int NT);
void gemm_5(int* A,int* B,long* C,int N,int NT);
void gemm_6(int* A,int* B,long* C,int N,int NT);
void gemm_7(int* A,int* B,long* C,int N,int NT);
void gemm_8(int* A,int* B,long* C,int N,int NT);
void main(int argc,char** argv)
{
    if(argc!=4 && argc!=3)
    {
        printf("shoule param: ./exe  [(int)para1:select func] [(int)para2: dim] [(int)para3 thread_nums]\n");
        return;
    }
    int func=1;
    int N=12;
    int NT=1;
    if(argc==3)
    {
        func = atoi(argv[1]);
        N = atoi(argv[2]);
    }
    else if(argc==4)
    {
        func = atoi(argv[1]);
        N = atoi(argv[2]);
        NT = atoi(argv[3]);
    }
 
    int r,c,k;
    long sum=0;
	int* A = calloc(N*N,sizeof(int));
	int* B = calloc(N*N,sizeof(int));
	long* C = calloc(N*N,sizeof(long));
    // 初始化
    for(int r=0;r<N;r++)
    {
        for(int c=0;c<N;c++)
            A[r*N+c] = B[r*N+c] = r+c;
    }
 
    switch (func)
    {
    case 1:
        NT = 1;
        gemm_1(A,B,C,N,NT);
        break;
    case 2:
        gemm_2(A,B,C,N,NT);
        break;
    case 3:
        gemm_3(A,B,C,N,NT);
        break;
    case 4:
        NT = 1;
        gemm_4(A,B,C,N,NT);
        break;
    case 5: //忽略，无法 collapse(5)。计算错误
        // gemm_5(A,B,C,N,NT);
        break;
    case 6:
        NT = 1;
        gemm_6(A,B,C,N,NT);
        break;
    case 7:
        gemm_7(A,B,C,N,NT);
        break;
    case 8:
        gemm_8(A,B,C,N,NT);
        break;
    default:
        break;
    }
 
	free(A);
	free(B);
	free(C);
}
 
long sum(long* C,int n)
{
	long sum=0;
	for(int i=0;i<n*n;i++)
		sum+=C[i];
	return sum;
}
 
void print_matrix(long* matrix, int rows, int cols)
{
	printf("matrix:\n");
	for (size_t i = 0; i < rows; i++)
	{
		for (size_t j = 0; j < cols; j++)
			printf("%ld\t",matrix[i*cols+j]);
		printf("\n");
	}
}
 
void print_matrixI(int* matrix, int rows, int cols)
{
	printf("matrix:\n");
	for (size_t i = 0; i < rows; i++)
	{
		for (size_t j = 0; j < cols; j++)
			printf("%d\t",matrix[i*cols+j]);
		printf("\n");
	}
}

(1) 矩阵不分块 串行

void gemm_1(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    for(r=0;r<N;r++)//A 行遍历
    {
        for(c=0;c<N;c++)//B 列遍历
        {
            // for循环 变量规约
            long sum=0;
            for(k=0;k<N;k++)//A B K方向遍历
                sum += A[r*N+k] * B[k*N+c];
            C[r*N+c] = sum;
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

(2) 矩阵不分块 并行 collapse(2)

void gemm_2(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    omp_set_num_threads(NT);  
    // 使用collapse(3) C[r][c] 的不同数据再不同进程之间，会造成数据竞争，需要原子操作,造成速度过慢。   
    // 使用collapse(3) for 默认循环指标变量私有，需要private(k)
    #pragma omp parallel for collapse(2) schedule(guided) proc_bind(close)
    for(int r=0;r<N;r++)//A 行遍历
    {
        for(int c=0;c<N;c++)//B 列遍历
        {
            #if 0   //临界区+collapse(3)
            for(int k=0;k<N;k++)//A B K方向遍历
            {
                #pragma omp critical
                    C[r*N+c] += A[r*N+k] * B[k*N+c];
            }
            #else
            long sum=0; //使用局部变量减少传入参数引用
            for(int k=0;k<N;k++)//A B K方向遍历
                sum += A[r*N+k] * B[k*N+c];
            C[r*N+c] = sum;
            #endif
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrixI(A,N,N);
    // printf("-----\n");
    // print_matrix(C,N,N);
 
}

(3) 矩阵不分块 并行 collapse(2)+ reduction

void gemm_3(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    omp_set_num_threads(NT);
    #pragma omp parallel for collapse(2) schedule(guided) proc_bind(close)
    for(r=0;r<N;r++)//A 行遍历
    {
        for(c=0;c<N;c++)//B 列遍历
        {
            // for循环 变量规约 默认k私有
            long sum=0;
            #pragma omp parallel for reduction(+:sum)  schedule(guided) proc_bind(close)
            for(k=0;k<N;k++)//A B K方向遍历
                sum += A[r*N+k] * B[k*N+c];
            C[r*N+c] = sum;
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

(4) 矩阵一维行列分块 A块 4N，B块 N4，串行，K在外层

void gemm_4(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    for(r=0;r<N;r+=4)//A 行遍历
    {
        for(c=0;c<N;c+=4)//B 列遍历
        {
            for(int k=0;k<N;k++) //K方向遍历
            {
                for(int nr=0;nr<4;nr++) //A块内行遍历
                {
                    for(int nc=0;nc<4;nc++) //B快内列遍历
                        C[(r+nr)*N+c+nc] += A[(r+nr)*N+k] * B[k*N+c+nc];
                }
            }
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

(5) 矩阵一维行列分块 A块 4N，B块 N4，并行，K在外层，collapse。2层外 3层内循环，无法 collapse(5)。

void gemm_5(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    omp_set_num_threads(NT);
    #pragma omp parallel for collapse(5) schedule(guided) proc_bind(close)
    for(r=0;r<N;r+=4)//A 行遍历
    {
        for(c=0;c<N;c+=4)//B 列遍历
        {
            for(int k=0;k<N;k++) //K方向遍历
            {
                for(int nr=0;nr<4;nr++) //A块内行遍历
                {
                    for(int nc=0;nc<4;nc++) //B快内列遍历
                    {
                        C[(r+nr)*N+c+nc] += A[(r+nr)*N+k] * B[k*N+c+nc];
                    }
                }
            }
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

(6) 矩阵一维行列分块 A块 4N，B块 N4，串行，K在内层。

void gemm_6(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    for(r=0;r<N;r+=4)//A 行遍历
    {
        for(c=0;c<N;c+=4)//B 列遍历
        {
            for(int nr=0;nr<4;nr++) //A块内行遍历
            {
                for(int nc=0;nc<4;nc++) //B快内列遍历
                {
                    long sum =0;
                    for(int k=0;k<N;k++) //K方向遍历
                        sum += A[(r+nr)*N+k] * B[k*N+c+nc];
                    C[(r+nr)*N+c+nc] = sum;
                }
            }
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

(7) 矩阵一维行列分块 A块 4N，B块 N4，并行，K在内层，collapse(2 + 3)

void gemm_7(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    omp_set_num_threads(NT);
    #pragma omp parallel for collapse(2) schedule(guided) proc_bind(close)
    for(r=0;r<N;r+=4)//A 行遍历
    {
        for(c=0;c<N;c+=4)//B 列遍历
        {
            #pragma omp parallel for collapse(3) schedule(guided) proc_bind(close)
            for(int nr=0;nr<4;nr++) //A块内行遍历
            {    
                for(int nc=0;nc<4;nc++) //B快内列遍历
                {
                    for(int k=0;k<N;k++) //K方向遍历
                        C[(r+nr)*N+c+nc] += A[(r+nr)*N+k] * B[k*N+c+nc];
                }
            }
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

(8) 矩阵一维行列分块 A块 4N，B块 N4，并行，K在内层，collapse (2+2)+ reduction

void gemm_8(int* A,int* B,long* C,int N,int NT)
{
    struct timeval start,end;
    float time;
    int r,c,k;
    //串行计算代码
    gettimeofday(&start,NULL);  //开始时间
    omp_set_num_threads(NT);
    #pragma omp parallel for collapse(2) private(r,c) schedule(guided) proc_bind(close)
    for(r=0;r<N;r+=4)//A 行遍历
    {
        for(c=0;c<N;c+=4)//B 列遍历
        {
            #pragma omp parallel for collapse(2) schedule(guided) proc_bind(close)
            for(int nr=0;nr<4;nr++) //A块内行遍历
            {
                for(int nc=0;nc<4;nc++) //B快内列遍历
                {
                    long sum =0;
                    #pragma omp parallel for reduction(+:sum)  schedule(guided) proc_bind(close)
                    for(int k=0;k<N;k++) //K方向遍历
                        sum += A[(r+nr)*N+k] * B[k*N+c+nc];
                    C[(r+nr)*N+c+nc] = sum;
                }
            }
        }
    }
    gettimeofday(&end,NULL);
    time = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
    printf("func %s N %d threads_num %d sum %ld useSec %f\n",__func__,N,NT,sum(C,N),time);
    // print_matrix(C,N,N);
}

4. 测试数据

串行耗时（秒） FIELD2 FIELD3 FIELD4 FIELD5
函数名 线程数\维度 128 512 1024
gemm_1 矩阵不分块 串行 1 0.01786 1.52600 47.10980
gemm_4 矩阵一维行列分块 A块 4N，B块 N4，串行，K在外层 1 0.03303 2.13628 22.47993
gemm-6 矩阵一维行列分块 A块 4N，B块 N4，串行，K在内层， 1 0.02078 1.64609 65.69423
并行耗时（秒）
函数名 线程数\维度 128 512 1024
gemm_2 矩阵不分块 并行 collapse(2) 4 0.00515 0.36207 13.16278
8 0.00303 0.18024 6.34872
16 0.00382 0.12891 3.39717
32 0.00471 0.10166 2.02269
gemm_3 矩阵不分块 并行 collapse(2) + reduction 4 0.01050 0.54038 11.69852
8 0.00593 0.23172 6.51897
16 0.00807 0.16333 3.60325
32 0.01048 0.12297 2.12352
gemm_7 矩阵一维行列分块 A块 4N，B块 N4，并行，K在内层，collapse(2+3) 4 0.00893 0.71235 18.41659
8 0.00509 0.33909 8.78864
16 0.00803 0.23421 4.56388
32 0.01073 0.15445 2.65277
gemm_8 矩阵一维行列分块 A块 4N，B块 N4，并行，K在内层，collapse(2+2) + reduction 4 0.01120 0.58372 16.70111
8 0.00622 0.26817 7.68263
16 0.00756 0.14903 4.24847
32 0.01052 0.13294 2.71985

5. 结果分析

(1). 关于线程数：较小的数据规模应用较低的线程数，否则线程创建销毁影响效率。较大的数据规模应采用更大的线程数量。总之数据规模应该与线程数量匹配才能效率较大。

(2). 关于reduction：效率没有开启向量化优化 O1及以上有效（测试了1024维度4/32线程加速明显），在问题规模较小（512维度及以下）即使没有采用向量化优化也比 reduction 效率高。

(3). 关于 collapse 和 reduction：测试了 collapse(2+3) 和 collapse(2+2) + reduction，在数据规模较小及线程数少（128维度和8线程及以下）全collapse效率高，除此之外随着数据规模和线程数增加，collapse 的任务规模数量级增加，效率不如 reduction。

posted @   安洛8  阅读(55)  评论(0)  编辑  收藏  举报
相关博文：
· 高性能计算-gemm串行计算优化(3)
· 高性能计算-openmp-多线程缓存一致性(9)
· OpenMP实验
· OpenMP学习 第四章 线程与OpenMP编程模型
· 借助 mperf 进行矩阵乘法极致优化
阅读排行：
· DeepSeek 开源周回顾「GitHub 热点速览」
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布：重大改进与新特性概览！
· AI与.NET技术实操系列（二）：开始使用ML.NET
· 单线程的Redis速度为什么快？
点击右上角即可分享
微信分享提示
AI FOR CODE 大赛