按行分配
思路和MPI基本类似,不过OpenMP是共享内存的,不必做分发和聚集,申请的矩阵空间就不必是完全连续的。
1 #include<stdio.h> 2 #include<omp.h> 3 #include<stdlib.h> 4 5 #define N 400 //规模(方针的阶数) 6 int i,j;//通用游标 7 double **mat=NULL;//矩阵对象 8 int *vec=NULL;//列向量对象 9 double *result=NULL;//结果向量对象 10 11 int main(int argc,char* argv[]) 12 { 13 //读入10进制表示的进程数 14 int thrdCnt=strtol(argv[1],NULL,10); 15 //矩阵一级挂载空间 16 mat=(double**)malloc(N*sizeof(double *)); 17 //存向量的空间 18 vec=(int *)malloc(N*sizeof(int)); 19 //存结果的空间 20 result=(double *)malloc(N*sizeof(double)); 21 22 //并行for:申请存矩阵的空间 23 //因为OpenMP不需要分发聚集等,所以不必做连续空间申请 24 # pragma omp parallel for num_threads(thrdCnt) 25 for(i=0;i<N;i++) 26 mat[i]=(double*)malloc(N*sizeof(double)); 27 28 //并行for:为矩阵和列向量赋值(实际做时是从文件读入) 29 # pragma omp parallel for num_threads(thrdCnt) private(j) 30 for(i=0;i<N;i++) 31 { 32 vec[i]=1; 33 for(j=0;j<N;j++) 34 mat[i][j]=j; 35 } 36 37 //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断 38 # pragma omp parallel for num_threads(thrdCnt) 39 for(i=0;i<N;i++) 40 result[i]=mat[i][0]*vec[0]; 41 42 //并行for:计算结果 43 # pragma omp parallel for num_threads(thrdCnt) private(j) 44 for(i=0;i<N;i++) 45 for(j=1;j<N;j++) 46 result[i]+=mat[i][j]*vec[j]; 47 48 //采样输出结果看一下 49 for(i=0;i<N;i+=N/11) 50 printf("%f\n",result[i]); 51 52 return 0; 53 }
输出
1 [lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_row.o omp_row.c 2 [lzh@hostlzh OpenMP]$ ./omp_row.o 7 3 79800.000000 4 79800.000000 5 79800.000000 6 79800.000000 7 79800.000000 8 79800.000000 9 79800.000000 10 79800.000000 11 79800.000000 12 79800.000000 13 79800.000000 14 79800.000000 15 [lzh@hostlzh OpenMP]$
按列分配
按列分配有很多细节要注意,如果不对result数组保护则可能会发生冲突,如果用critical或者atomic会导致计算过程实际是串行的,虽然atomic的加解锁过程是critical的7倍。
1 #include<stdio.h> 2 #include<omp.h> 3 #include<stdlib.h> 4 5 #define N 400 //规模(方针的阶数) 6 int i,j;//通用游标 7 double **mat=NULL;//矩阵对象 8 int *vec=NULL;//列向量对象 9 double *result=NULL;//结果向量对象 10 11 int main(int argc,char* argv[]) 12 { 13 //读入10进制表示的进程数 14 int thrdCnt=strtol(argv[1],NULL,10); 15 //矩阵一级挂载空间 16 mat=(double**)malloc(N*sizeof(double *)); 17 //存向量的空间 18 vec=(int *)malloc(N*sizeof(int)); 19 //存结果的空间 20 result=(double *)malloc(N*sizeof(double)); 21 22 //并行for:申请存矩阵的空间 23 //因为OpenMP不需要分发聚集等,所以不必做连续空间申请 24 # pragma omp parallel for num_threads(thrdCnt) 25 for(i=0;i<N;i++) 26 mat[i]=(double*)malloc(N*sizeof(double)); 27 28 //并行for:为矩阵和列向量赋值(实际做时是从文件读入) 29 # pragma omp parallel for num_threads(thrdCnt) private(j) 30 for(i=0;i<N;i++) 31 { 32 vec[i]=1; 33 for(j=0;j<N;j++) 34 mat[i][j]=j; 35 } 36 37 //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断 38 # pragma omp parallel for num_threads(thrdCnt) 39 for(i=0;i<N;i++) 40 result[i]=mat[i][0]*vec[0]; 41 42 //创建存放每个线程临时结果的数组,初始化为0 43 //这里calloc还能像前面那样改进,但为了程序可读性暂时不做改进 44 double* sum; 45 46 //OpenMP默认块划分给每个进程(除了最后一个进程)的循环次数 47 //是总的循环次数处以进程数向上取整,需要先计算出这个数 48 int needPlus;//记录是否向上+1 49 needPlus=((N-1)%thrdCnt==0)?0:1; 50 //每个进程(除了最后一个进程)的循环次数 51 int numThrdFor=(N-1)/thrdCnt+needPlus; 52 53 //并行for:计算结果,按列分配时这个大的外层for用j,且从1开始 54 # pragma omp parallel for num_threads(thrdCnt) \ 55 private(i) firstprivate(sum)//对i只要每个线程私有,对sum数组还需初值 56 for(j=1;j<N;j++) 57 { 58 //本线程第一次运行时创建sum空间 59 if((j-1)%numThrdFor==0) 60 { 61 //创建存放自己线程临时结果的数组,初始化为0 62 //这里calloc还能像前面那样改进,但为了程序可读性暂时不做改进 63 sum=calloc(sizeof(double),N); 64 } 65 66 //对于这其中的每一短行 67 //(这个i被private保护,为每个线程单独创建了私有的i) 68 for(i=0;i<N;i++) 69 { 70 //如果用critical保护所有result[i]本质上这个计算完全是串行的 71 //想使用reduction子句,但是预编译只会做一次,而不会随for循环 72 //使用atomic使得这条语句是原子的,要执行就执行完而不拆分 73 //但用atomic这个计算也仍然是串行的,只是加解锁比critical快 74 //因此对每个线程拷贝分配一个求和数组,才能实现并行 75 sum[i]+=mat[i][j]*vec[j];//加到自己线程的sum数组上 76 } 77 78 //仅当本线程即将结束前,把算好的sum数组加上来 79 //判断条件:从1开始计数下能整除次数,或当前是最后一次循环 80 if(j%numThrdFor==0 || j==N-1) 81 { 82 for(i=0;i<N;i++) 83 # pragma omp atomic//在这里使用atomic保护result[i] 84 result[i]+=sum[i]; 85 } 86 } 87 88 //采样输出结果看一下 89 for(i=0;i<N;i+=N/11) 90 printf("%f\n",result[i]); 91 92 return 0; 93 }
输出
1 [lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_col.o omp_col.c 2 [lzh@hostlzh OpenMP]$ ./omp_col.o 12 3 79800.000000 4 79800.000000 5 79800.000000 6 79800.000000 7 79800.000000 8 79800.000000 9 79800.000000 10 79800.000000 11 79800.000000 12 79800.000000 13 79800.000000 14 79800.000000 15 [lzh@hostlzh OpenMP]$
按块分配
每次循环都会重新开个线程,虽然操作了共享变量但是测试时没出现问题(为啥)。
1 #include<stdio.h> 2 #include<omp.h> 3 #include<stdlib.h> 4 5 #define N 400 //规模(方针的阶数) 6 int i,j;//通用游标 7 double **mat=NULL;//矩阵对象 8 int *vec=NULL;//列向量对象 9 double *result=NULL;//结果向量对象 10 int blkSqrt=-1;//块数的开算数平方 11 12 //自己的求平方根的函数,因为math里的sqrt报错 13 int mysqrt(int k) 14 { 15 int i; 16 for(i=0;i*i<k;i++) 17 ; 18 if(i*i==k) 19 return i; 20 return -1; 21 } 22 23 int main(int argc,char* argv[]) 24 { 25 //读入10进制表示的进程数 26 int thrdCnt=strtol(argv[1],NULL,10); 27 //判断块数合法性 28 if(blkSqrt=mysqrt(thrdCnt)==-1)//块数不是完全平方数 29 { 30 printf("块数不是完全平方数!\n"); 31 return 0; 32 } 33 //矩阵一级挂载空间 34 mat=(double**)malloc(N*sizeof(double *)); 35 //存向量的空间 36 vec=(int *)malloc(N*sizeof(int)); 37 //存结果的空间 38 result=(double *)malloc(N*sizeof(double)); 39 40 //并行for:申请存矩阵的空间 41 //因为OpenMP不需要分发聚集等,所以不必做连续空间申请 42 # pragma omp parallel for num_threads(thrdCnt) 43 for(i=0;i<N;i++) 44 mat[i]=(double*)malloc(N*sizeof(double)); 45 46 //并行for:为矩阵和列向量赋值(实际做时是从文件读入) 47 # pragma omp parallel for num_threads(thrdCnt) private(j) 48 for(i=0;i<N;i++) 49 { 50 vec[i]=1; 51 for(j=0;j<N;j++) 52 mat[i][j]=j; 53 } 54 55 //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断 56 # pragma omp parallel for num_threads(thrdCnt) 57 for(i=0;i<N;i++) 58 result[i]=mat[i][0]*vec[0]; 59 60 //并行for:按块分配计算结果,即行列分别分开 61 # pragma omp parallel for num_threads(blkSqrt) 62 for(i=0;i<N;i++) 63 # pragma omp parallel for num_threads(blkSqrt) 64 for(j=1;j<N;j++) 65 result[i]+=mat[i][j]*vec[j]; 66 67 //采样输出结果看一下 68 for(i=0;i<N;i+=N/11) 69 printf("%f\n",result[i]); 70 71 return 0; 72 }
输出
1 [lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_blk.o omp_blk.c 2 [lzh@hostlzh OpenMP]$ ./omp_blk.o 9 3 79800.000000 4 79800.000000 5 79800.000000 6 79800.000000 7 79800.000000 8 79800.000000 9 79800.000000 10 79800.000000 11 79800.000000 12 79800.000000 13 79800.000000 14 79800.000000 15 [lzh@hostlzh OpenMP]$
标签:int,double,omp,79800.000000,OpenMP,NULL,方阵,乘法 From: https://www.cnblogs.com/ybqjymy/p/17462930.html