OpenMP(三)梯形积分法以及pi值估计

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

1.梯形积分法求解函数下方面积

求解函数下方的面积，可以将面积分成很多小的梯形条，计算每个梯形条，然后使用reduction求和。

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

/*
    利用梯形积分法计算函数下方的面积
    将面积分为很多小条梯形，然后求和
*/
#define f(x)  1/(1+x*x)     //首先随意定义一个函数，计算函数下方的面积
int main() {
    int thread_num = 10;    //线程数
    int n = 100000;         //梯形数
    double a = 0.0;         //积分下限
    double b = 100.0;       //积分上限
    double h = (b-a)/n;     //梯形宽度
    double sum = 0.0;       //梯形面积和

    //计算时间
    double start_time = omp_get_wtime();

    //并行计算
    #pragma omp parallel for reduction(+:sum) num_threads(thread_num)
    for(int i=0; i<n; i++) {
        sum += (f(a+i*h)+f(a+(i+1)*h))*h/2;
    }

    //计算时间
    double end_time = omp_get_wtime();
    double compute_time = end_time - start_time;

    //输出结果
    printf("sum = %lf\n", sum);
    printf("compute time = %lf\n", compute_time);
}

输出结果

sum = 4.615121
compute time = 0.001000

2.蒙特卡洛方法求解Pi的估计值

蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法。在正方形内切一个圆，那么圆的面积和正方形面积之比是pi/4，我们可以随机往正方形内掷色子，只要色子扔的够多，那么它的圆形内部的概率也刚好是pi/4，进而可以通过统计色子在圆内的概率来计算pi的近似值。

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

/*
    蒙特卡洛方法计算圆周率pi
*/
int main() {
    int thread_num = 20;
    int n = 1000000000;         //掷色子的数目
    int count = 0;          //落在圆内的点数
    double x, y;            //随机点坐标
    double rand_r = 1;      //半径
    unsigned seed = time(NULL);

    //计算时间
    double start_time = omp_get_wtime();

    #pragma omp parallel for reduction(+:count) num_threads(thread_num) private(x, y)
    for(int i=0; i<n; i+=1) {
        x = rand() / (double) RAND_MAX;
        y = rand() / (double) RAND_MAX;
        if(x*x+y*y<1)   count++;
    }

    double pi = 4 * (count / (double)n);

    double end_time = omp_get_wtime();
    double compute_time = end_time - start_time;

    printf("pi = %lf\n", pi);
    printf("compute time = %lf\n", compute_time);
    return 0;
}

运行结果

pi = 3.141356
compute time = 6.462000

如果将thread_num从20改为10，那么运行结果将为

pi = 3.141353
compute time = 7.477000

每次的计算时间都不同，但基本上是比线程数为20所花时间长

但如果将thread_num从20提升为32，运行结果如下

pi = 3.141399
compute time = 6.653000

可以看出，并不是线程数越大越好，我们可以按照2×核心数-1设置。代码部分对应修改如下：

int thread_num = 2*omp_get_num_procs()-1;

运行结果：

pi = 3.141343
compute time = 6.647000

好吧，并没有比线程数为20好到哪里去(捂脸)。如果实在要找到最优线程数，可以循环计算不同线程数下的计算时间。