OpenMP并行计算PI

程序包括多个函数：

1. 单线程计算PI
2. 开启并行域并行（SPMD并行模式）
3. 使用for制导指令
4. 使用reduction子句
5. 使用critical子句

---

#include <stdio.h>
#include <time.h>
#include <omp.h>
#include <math.h>
#define realPI acos(-1)
const int num_steps = 300000000;
const int NUM_THREADS = 20;

//单线程计算Pi
void single()
{
    puts("single thread...");
    double step, x, pi, sum = 0.0, timetot;
    step = 1.0 / (double)num_steps;
    clock_t start, end;
    start = clock();
    int i = 1;
    for (; i <= num_steps; i++)
    {
        x = (i + 0.5) * step;
        sum += 4.0 / (1.0 + x * x);
    }
    pi = step * sum;
    end = clock();
    timetot = (double)(end - start) / CLOCKS_PER_SEC;
    printf("Pi=%.20f  Running time=%f s\n", pi, timetot);
    puts("");
}

// 并行域并行（SPMD并行模式）
// void SPMD()
// {
//     puts("SPMD mode...");
//     int i;
//     double step, start_time, end_time, pi = 0.0, sum[NUM_THREADS];
//     step = 1.0 / (double)num_steps;
//     omp_set_num_threads(NUM_THREADS);
//     start_time = omp_get_wtime();
// #pragma omp parallel
//     {
//         double x;
//         int id = omp_get_thread_num();
//         for (i = id, sum[id] = 0.0; i < num_steps; i = i + NUM_THREADS)
//         {
//             x = (i + 0.5) * step;
//             sum[id] += 4.0 / (1.0 + x * x);
//         }
//     }
//     for (i = 0; i < NUM_THREADS; i++)
//         pi += sum[i] * step;
//     end_time = omp_get_wtime();
//     printf("Pi=%.20f  Running time=%f s\n", pi, end_time - start_time);
//     puts("");
// }

// 使用for制导指令
void withfor()
{
    puts("using '#pragma omp for'...");
    int i;
    double step, start_time, end_time, pi = 0.0, sum[NUM_THREADS];
    step = 1.0 / (double)num_steps;
    omp_set_num_threads(NUM_THREADS);
    start_time = omp_get_wtime();
#pragma omp parallel private(i)
    {
        int id = omp_get_thread_num();
        sum[id] = 0.0;
        double x;
#pragma omp for
        for (i = 0; i < num_steps; i++)
        {
            x = (i + 0.5) * step;
            sum[id] += 4.0 / (1.0 + x * x);
        }
    }
    for (i = 0; i < NUM_THREADS; i++)
        pi += sum[i] * step;
    end_time = omp_get_wtime();
    printf("Pi=%.20f  Running time=%f s\n", pi, end_time - start_time);
    puts("");
}

// 使用reduction子句
void reduction()
{
    puts("using '#pragma omp parallel for reduction(+:sum)'...");
    int i;
    double x, sum = 0.0, step, start_time, end_time, pi = 0.0;
    step = 1.0 / (double)num_steps;
    omp_set_num_threads(NUM_THREADS);
    start_time = omp_get_wtime();
#pragma omp parallel for reduction(+ \
                                   : sum) private(x)
    for (i = 0; i < num_steps; i++)
    {
        x = (i + 0.5) * step;
        sum = sum + 4.0 / (1.0 + x * x);
    }
    pi = step * sum;
    end_time = omp_get_wtime();
    printf("Pi=%.20f  Running time=%f s\n", pi, end_time - start_time);
    puts("");
}

// 使用critical子句
void critical()
{
    puts("using '#pragma omp critical'...");
    int i;
    double x, sum = 0.0, step, start_time, end_time, pi = 0.0;
    step = 1.0 / (double)num_steps;
    omp_set_num_threads(NUM_THREADS);
    start_time = omp_get_wtime();
#pragma omp parallel private(x, i, sum)
    {
        int id = omp_get_thread_num();
        for (i = id, sum = 0.0; i < num_steps; i += NUM_THREADS)
        {
            x = (i + 0.5) * step;
            sum += 4.0 / (1.0 + x * x);
        }
//critical指定某一区域的代码，每次只能同时被一个线程执行。
#pragma omp critical
        pi += sum * step;
    }
    end_time = omp_get_wtime();
    printf("Pi=%.20f  Running time=%f s\n", pi, end_time - start_time);
    puts("");
}

int main()
{
    printf("Caculating PI(%.20f) in 300000000 iterations with %d threads...\n", realPI, NUM_THREADS);
    puts("");
    single();
    // SPMD();
    withfor();
    reduction();
    critical();
    return 0;
}