OpenMP入门:求pi

思路:应用积分的思想求解4/(x2+1)在0-1区间积分的值。将0-1积分区间分成n份,使n尽可能大。这时,被划分的每个小区域就可近似为一个小矩形,其宽为1/n,长可取小区域内的任意高度。求出每个小矩形的面积相加后即得到π值。

在程序中,共n=4个线程,第i个线程负责计算[i/n,(i+1)/n]部分,各线程维护自己的部分和,最后求和得到总的部分和。

OpenMP 是非常方便的

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

#define int64_t long long

const int64_t nThreads = 4;
const int64_t num_steps = 1e9;
double step;
double answer_pi, answer = 0.0;
double startTime, runTime;
double sum[nThreads];

int main()
{
    step = 1.0 / (double)num_steps;

    omp_set_num_threads(nThreads);
    answer = 0.0;
    startTime = omp_get_wtime();

#pragma omp parallel
    {
        int64_t id = omp_get_thread_num();
        int64_t numthreads = omp_get_num_threads();
        double x;
        double sum = 0;
        for (int64_t i = id * num_steps / nThreads; i < (id + 1) * num_steps / nThreads; i++)
        {
            x = (i + 0.5) * step;
            sum += 4.0 / (1.0 + x * x);
        }
#pragma omp critical
        answer += sum;
    }

    answer_pi = step * answer;
    runTime = omp_get_wtime() - startTime;
    printf("pi=%.14lf
time=%.10lf sec
%lld threads
", answer_pi, runTime, nThreads);
}