OpenMP-筛选法求素数

开启掘金成长之旅！这是我参与「掘金日新计划 · 12 月更文挑战」的第20天，点击查看活动详情

写在前面

上一篇博客使用MPI实现筛选法求素数，这只是并行计算一个简单案例，当计算量比较大(例如求10^100以内的质数)时，MPI的优势就体现出来了：先将任务分发出去，然后每个处理器并行计算自己区域的质数，最后通过MPI收集每个处理器计算的结果。其实，在单机内也有加速的方法，通过多个线程并发处理数据 (注意是并发)。，这就是OpenMP，它通过指导语句能方便地启动多个线程来进行计算，对于循环计算尤其友好。

回顾朴素的循环方法

#include<stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
void main()
{
    int N = 100;
    char* identifies = (char*)malloc(sizeof(char)*N);
    memset(identifies, 0, N);
    int i,j,k;
    for(i=2;i<=(int)(sqrt(N));i++)     //从2-10作为最小值开始筛选
    {
        for(j=2;j*i<N;j++)   //找出最小数的倍数
        {
            identifies[i*j]=1;     //把最小数的倍数赋值1
        }
    }
    for(k=2;k<N;k++)   //因为1不是素数，所以从2开始遍历输出素数
    {
        if(identifies[k]==0)       //在a数组中如果数组元素是0，那么相应下标就是素数
        {
            printf("%4d",k);
        }
    }
}

我们只需要在计算部分加上制导语句即可：

sushu.c

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

int main(int argc, char** argv) {
    int local_size, index;
    int n_numbers = 32;

    int nthreads = 4;
    omp_set_num_threads(nthreads);
    int tid = omp_get_num_threads();

    char* mark = (char*)malloc(sizeof(char)*n_numbers);         //使用char节省空间
    memset(mark, 0, n_numbers);

    #pragma omp parallel for
    for(int k=2; k<=(int)(pow((double)(n_numbers), 0.5)); k++) {
        int max = n_numbers/k;
        for(int j=2; j<=max; j++) {
            mark[j*k-1] = 1;
            printf("%d ", j*k);
        }
    }
    
    int count = 0;
    printf("\nprime numbers are: \n");
    for (int i = 0; i < n_numbers; i++) {
        count += 1- mark[i];
        if(!mark[i]&& (i)!=0)    printf("%d ", i+1);
    }
    printf("\n");
    return 0;
}

注意几点：制导语句下面的for嵌套循环中需要修改：

原：

#pragma omp parallel for collapse(2) private(k) shared(mark, n_numbers)
    for(k=2; k<=(int)(pow((double)(n_numbers), 0.5)); k++) {
        for(j=2; j*k<=n_numbers; j++) {
            mark[j*k-1] = 1;
            printf("%d ", j*k);
        }
    }

这段代码使用了循环展开，实际上该代码会报错：

error: invalid controlling predicate

主要原因是：for(j=2; j*k<=n_numbers; j++)中的j和k不能同时出现。故，只能放弃collapse，改用sushu.c的代码。

运行结果：

6 9 12 15 18 21 24 27 30 10 15 20 25 30 8 12 16 20 24 28 32 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 
prime numbers are:
2 3 5 7 11 13 17 19 23 29 31

代码中设置的线程数目为4，它会自动为每个线程分配计算空间，也就是每个线程处理25个数。当计算量到一定程度，不同的线程能够并发操作分担任务，提升计算速度。