用OpenMP对矩阵乘法运算进行优化

用OpenMP对矩阵乘法运算进行优化

数学算式：


算法流程伪代码：

在进行优化之前，先来讨论一下
用parallel for处理循环易遇到问题：
参考文章：openmp中parallel for处理双重循环遇到问题的分析和解决方案
作为初学者使用parallel for的双重循环，代码结果一直计算不对，但是一直没找到原因。后来发现原因挺简单的，就是在双重循环的时候，假如第一层是遍历 i，第二层是遍历 j，这个 j 如果是omp_set_num_threads()之前声明的，就说明 j 被多个线程共用，那么可能这个线程正在使用 j，然后被另一个线程更改了
例如：

#include <omp.h>
#include <stdio.h>
#include <iostream>
using namespace std;
int main()
{
    int i, j;
    printf("initila i address is %p, initial j address is %p\n", &i, &j);
    
    omp_set_num_threads(2);

    #pragma omp parallel for
    for (i = 0; i < 4; ++i)
    {
        for (j = 0; j < 1; ++j)
        {
            printf("thread is %d, i=%d, i address is %p, j=%d, j address is %p\n", 
                    omp_get_thread_num(), i, &i, j, &j);
        }
    }

    return 0;
}

initila i address is 000000000070FE0C, initial j address is 000000000070FE08
thread is 0, i=0, i address is 000000000070FD4C, j=0, j address is 000000000070FE08
thread is 0, i=1, i address is 000000000070FD4C, j=0, j address is 000000000070FE08
thread is 1, i=2, i address is 000000000141FDEC, j=0, j address is 000000000070FE08
thread is 1, i=3, i address is 000000000141FDEC, j=0, j address is 000000000070FE08

可以看到，一共两个线程，i 和 j 都是在开辟线程之前定义，但是openmp parallel for会为每个线程新创建一个 i，但是并不会新建一个 j，j 还是最开始创建的那个。这就导致双重循环中，第二重循环的 j 产生数据冲突

解决方案：
使用parallel for，只会把最外层的循环分给不同的线程，同一个外层循环的内层的循环只能由一个线程执行，要想实现同一个外层循环的内层循环也被多个线程执行，我们可以写循环的时候，内层循环写成：for (int j = 0; j < xxx; ++j)，就可以为每个线程创建不同的 j 了，避免数据冲突。或者我们可以把#pragma omp parallel for(private j)，强制让每个线程新创建 j

——————————————————
这里采用前者进行编写。

#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
#define N 300
int main()
{
	//随机生成300阶矩阵
	int a[N][N],b[N][N],c[N][N]={0};
	double t0,t1;
	for(int i=0;i<N;i++)
	{
		for(int j=0;j<N;j++)
		{
			a[i][j]=rand()%100;
			b[i][j]=rand()%100;
		}
	}
	
    //初始矩阵乘法运算
	t0=omp_get_wtime();
	for(int i=0;i<N;i++)
		for(int j=0;j<N;j++)
			for(int k=0;k<N;k++)
				c[i][j]+=a[i][k]*b[k][j];
	t1=omp_get_wtime();
	printf("初始矩阵乘法运行耗时：%f ms\n",(t1-t0)*1000); 
	
	//优化矩阵乘法运算 
	t0=omp_get_wtime();
	#pragma omp parallel for
	for(int i=0;i<N;i++)
		for(int j=0;j<N;j++)
            for(int k=0;k<N;k++)
				c[i][j]+=a[i][k]*b[k][j];
	t1=omp_get_wtime();
	printf("优化矩阵乘法（默认调度）运行耗时：%f ms\n",(t1-t0)*1000); 
	
	//优化矩阵乘法运算 （静态调度） 
	t0=omp_get_wtime();
	#pragma omp parallel for schedule(static)
	for(int i=0;i<N;i++)
		for(int j=0;j<N;j++)
            for(int k=0;k<N;k++)
				c[i][j]+=a[i][k]*b[k][j];
	t1=omp_get_wtime();
	printf("优化矩阵乘法（静态调度）运行耗时：%f ms\n",(t1-t0)*1000); 
	
	// 优化矩阵乘法运算 （动态调度） 
	t0=omp_get_wtime();
	#pragma omp parallel for schedule(dynamic)
	for(int i=0;i<N;i++)
		for(int j=0;j<N;j++)
            for(int k=0;k<N;k++)
				c[i][j]+=a[i][k]*b[k][j];
	t1=omp_get_wtime();
	printf("优化矩阵乘法（动态调度）运行耗时：%f ms\n",(t1-t0)*1000); 
	return 0;
}

结果：

这里可以看到static和dynamic差距并不大，应该是矩阵乘法任务本身for循环中的计算量大致相同的原因

对于稀疏矩阵，可以采用CSR法对数据进行储存

算法方面呢，也可以采用分块矩阵，提高运算效率：矩阵乘法的分块优化
分块的矩阵乘法是为了利用局部性提高cache（高速缓存）利用率

t0=omp_get_wtime();
	const int B=10;
	for(int i1=0;i1<N;i1+=B)
	{
		for(int j1=0;j1<N;j1+=B)
		{
			for(int k1=0;k1<N;k1+=B)
			{
				for(int i2=i1;i2<i1+B;i2++)
				{
					for(int j2=j2;j2<j1+B;j2++)
					{
						for(int k2=k1;k2<k1+B;k2++)
						c[i2][j2]+=a[i2][k2]*b[k2][j2]; 
					}
						
				}
					
			}
				
		}
			
	}
	t1=omp_get_wtime();
	printf("分块矩阵乘法运行耗时：%f ms\n",(t1-t0)*1000); 

将之与普通矩阵乘法进行对比，结果如下