【OpenMP学习笔记】4：传统形式的方阵向量并行乘法

按行分配

思路和MPI基本类似，不过OpenMP是共享内存的，不必做分发和聚集，申请的矩阵空间就不必是完全连续的。

#include<stdio.h>
#include<omp.h>
#include<stdlib.h>

#define N 400 //规模(方针的阶数)
int i,j;//通用游标
double **mat=NULL;//矩阵对象
int *vec=NULL;//列向量对象
double *result=NULL;//结果向量对象

int main(int argc,char* argv[])
{
    //读入10进制表示的进程数
    int thrdCnt=strtol(argv[1],NULL,10);
    //矩阵一级挂载空间
    mat=(double**)malloc(N*sizeof(double *));
    //存向量的空间
    vec=(int *)malloc(N*sizeof(int));
    //存结果的空间
    result=(double *)malloc(N*sizeof(double));

    //并行for:申请存矩阵的空间
    //因为OpenMP不需要分发聚集等,所以不必做连续空间申请
#   pragma omp parallel for num_threads(thrdCnt)
    for(i=0;i<N;i++)
        mat[i]=(double*)malloc(N*sizeof(double));   

    //并行for:为矩阵和列向量赋值(实际做时是从文件读入)
#   pragma omp parallel for num_threads(thrdCnt) private(j)
    for(i=0;i<N;i++)
    {
        vec[i]=1;
        for(j=0;j<N;j++)
            mat[i][j]=j;
    }

    //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断
#   pragma omp parallel for num_threads(thrdCnt)
    for(i=0;i<N;i++)
        result[i]=mat[i][0]*vec[0];

    //并行for:计算结果
#   pragma omp parallel for num_threads(thrdCnt) private(j)
    for(i=0;i<N;i++)
        for(j=1;j<N;j++)
            result[i]+=mat[i][j]*vec[j];    

    //采样输出结果看一下
    for(i=0;i<N;i+=N/11)
        printf("%f\n",result[i]);

    return 0;
}

输出

[lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_row.o omp_row.c
[lzh@hostlzh OpenMP]$ ./omp_row.o 7
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
[lzh@hostlzh OpenMP]$

按列分配

按列分配有很多细节要注意，如果不对result数组保护则可能会发生冲突，如果用critical或者atomic会导致计算过程实际是串行的，虽然atomic的加解锁过程是critical的7倍。

#include<stdio.h>
#include<omp.h>
#include<stdlib.h>

#define N 400 //规模(方针的阶数)
int i,j;//通用游标
double **mat=NULL;//矩阵对象
int *vec=NULL;//列向量对象
double *result=NULL;//结果向量对象

int main(int argc,char* argv[])
{
    //读入10进制表示的进程数
    int thrdCnt=strtol(argv[1],NULL,10);
    //矩阵一级挂载空间
    mat=(double**)malloc(N*sizeof(double *));
    //存向量的空间
    vec=(int *)malloc(N*sizeof(int));
    //存结果的空间
    result=(double *)malloc(N*sizeof(double));

    //并行for:申请存矩阵的空间
    //因为OpenMP不需要分发聚集等,所以不必做连续空间申请
#   pragma omp parallel for num_threads(thrdCnt)
    for(i=0;i<N;i++)
        mat[i]=(double*)malloc(N*sizeof(double));   

    //并行for:为矩阵和列向量赋值(实际做时是从文件读入)
#   pragma omp parallel for num_threads(thrdCnt) private(j)
    for(i=0;i<N;i++)
    {
        vec[i]=1;
        for(j=0;j<N;j++)
            mat[i][j]=j;
    }

    //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断
#   pragma omp parallel for num_threads(thrdCnt)
    for(i=0;i<N;i++)
        result[i]=mat[i][0]*vec[0];

    //创建存放每个线程临时结果的数组,初始化为0
    //这里calloc还能像前面那样改进,但为了程序可读性暂时不做改进
    double* sum;

    //OpenMP默认块划分给每个进程(除了最后一个进程)的循环次数
    //是总的循环次数处以进程数向上取整,需要先计算出这个数
    int needPlus;//记录是否向上+1
    needPlus=((N-1)%thrdCnt==0)?0:1;
    //每个进程(除了最后一个进程)的循环次数
    int numThrdFor=(N-1)/thrdCnt+needPlus;

    //并行for:计算结果,按列分配时这个大的外层for用j,且从1开始
#   pragma omp parallel for num_threads(thrdCnt) \
    private(i) firstprivate(sum)//对i只要每个线程私有,对sum数组还需初值
    for(j=1;j<N;j++)
    {
        //本线程第一次运行时创建sum空间
        if((j-1)%numThrdFor==0)
        {
            //创建存放自己线程临时结果的数组,初始化为0
            //这里calloc还能像前面那样改进,但为了程序可读性暂时不做改进
            sum=calloc(sizeof(double),N);
        }

        //对于这其中的每一短行
        //(这个i被private保护,为每个线程单独创建了私有的i)
        for(i=0;i<N;i++)
        {
            //如果用critical保护所有result[i]本质上这个计算完全是串行的
            //想使用reduction子句,但是预编译只会做一次,而不会随for循环
            //使用atomic使得这条语句是原子的,要执行就执行完而不拆分
            //但用atomic这个计算也仍然是串行的,只是加解锁比critical快
            //因此对每个线程拷贝分配一个求和数组,才能实现并行
            sum[i]+=mat[i][j]*vec[j];//加到自己线程的sum数组上
        }

        //仅当本线程即将结束前,把算好的sum数组加上来
        //判断条件:从1开始计数下能整除次数,或当前是最后一次循环
        if(j%numThrdFor==0 || j==N-1)
        {
            for(i=0;i<N;i++)
#               pragma omp atomic//在这里使用atomic保护result[i]
                result[i]+=sum[i];
        }
    }

    //采样输出结果看一下
    for(i=0;i<N;i+=N/11)
        printf("%f\n",result[i]);

    return 0;
}

输出

[lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_col.o omp_col.c
[lzh@hostlzh OpenMP]$ ./omp_col.o 12
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
[lzh@hostlzh OpenMP]$

按块分配

每次循环都会重新开个线程，虽然操作了共享变量但是测试时没出现问题(为啥)。

#include<stdio.h>
#include<omp.h>
#include<stdlib.h>

#define N 400 //规模(方针的阶数)
int i,j;//通用游标
double **mat=NULL;//矩阵对象
int *vec=NULL;//列向量对象
double *result=NULL;//结果向量对象
int blkSqrt=-1;//块数的开算数平方

//自己的求平方根的函数,因为math里的sqrt报错
int mysqrt(int k)
{
    int i;
    for(i=0;i*i<k;i++)
        ;
    if(i*i==k)
        return i;
    return -1;
}

int main(int argc,char* argv[])
{
    //读入10进制表示的进程数
    int thrdCnt=strtol(argv[1],NULL,10);
    //判断块数合法性
    if(blkSqrt=mysqrt(thrdCnt)==-1)//块数不是完全平方数
    {
        printf("块数不是完全平方数!\n");
        return 0;
    }
    //矩阵一级挂载空间
    mat=(double**)malloc(N*sizeof(double *));
    //存向量的空间
    vec=(int *)malloc(N*sizeof(int));
    //存结果的空间
    result=(double *)malloc(N*sizeof(double));

    //并行for:申请存矩阵的空间
    //因为OpenMP不需要分发聚集等,所以不必做连续空间申请
#   pragma omp parallel for num_threads(thrdCnt)
    for(i=0;i<N;i++)
        mat[i]=(double*)malloc(N*sizeof(double));   

    //并行for:为矩阵和列向量赋值(实际做时是从文件读入)
#   pragma omp parallel for num_threads(thrdCnt) private(j)
    for(i=0;i<N;i++)
    {
        vec[i]=1;
        for(j=0;j<N;j++)
            mat[i][j]=j;
    }

    //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断
#   pragma omp parallel for num_threads(thrdCnt)
    for(i=0;i<N;i++)
        result[i]=mat[i][0]*vec[0];

    //并行for:按块分配计算结果,即行列分别分开
#   pragma omp parallel for num_threads(blkSqrt)
    for(i=0;i<N;i++)
#       pragma omp parallel for num_threads(blkSqrt)
        for(j=1;j<N;j++)
            result[i]+=mat[i][j]*vec[j];

    //采样输出结果看一下
    for(i=0;i<N;i+=N/11)
        printf("%f\n",result[i]);

    return 0;
}

输出

[lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_blk.o omp_blk.c
[lzh@hostlzh OpenMP]$ ./omp_blk.o 9
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
79800.000000
[lzh@hostlzh OpenMP]$