使用Pthreads实现梯形法则求面积

问题

编写Pthreads程序实现梯形法则求面积

• 使用共享变量对线程的计算结果进行累加
• 使用busy-waiting，mutexes和semaphores实现对临界区域的互斥
• 根据结果分析每种方法的优缺点

思路

梯形法则求面积的数学知识：


梯形法则求面积的并行程序的伪代码：

主要代码（详细代码请点击查看附件）

1. busy-waiting
   我们使用了以下的全局变量：

long thread_count;//in
long long n;//in
double a;//left_endpt  in
double b;//right_endpt  in
long double h;
long long trap_count;

long double sum;
int flag;

在主线程中对flag初始化为0，使用flag来忙等待，控制进入临界区。以下是主函数的实现。

int main(int argc, char* argv[])
{
    long       thread;  /* Use long in case of a 64-bit system */
    pthread_t* thread_handles;
    double start, finish, elapsed;

    /* Get number of threads from command line */
    Get_args(argc, argv);

    thread_handles = (pthread_t*) malloc (thread_count*sizeof(pthread_t));
    sum = 0.0;
    flag = 0;
    h = (long double)(b-a)/n;
    trap_count = n/thread_count;

    GET_TIME(start);
    for (thread = 0; thread < thread_count; thread++)
        pthread_create(&thread_handles[thread], NULL, Thread_sum, (void*)thread);

    for (thread = 0; thread < thread_count; thread++)
        pthread_join(thread_handles[thread], NULL);
    GET_TIME(finish);
    elapsed = finish - start;

    sum = h*sum;
    printf("With n = %lld terms,\n", n);
    printf("   Our estimate of area = %.15Lf\n", sum);
    printf("The elapsed time is %e seconds\n", elapsed);

    GET_TIME(start);
    sum = Serial_trap(n);
    GET_TIME(finish);
    elapsed = finish - start;
    printf("   Single thread est    = %.15Lf\n", sum);
    printf("The elapsed time is %e seconds\n", elapsed);
    printf("          accurate area = 2.666666666666667\n");

    free(thread_handles);
    return 0;
}  /* main */

以下是用到的主要函数：

/*To calculate the trapezoidal area by multi-thread*/
void* Thread_sum(void* rank);

/* Only executed by main thread */
void Get_args(int argc, char* argv[]);
void Usage(char* prog_name);
double Serial_trap(long long n);

/* Function we're integrating */
double f(double x);

以下是线程函数的实现：

/*------------------------------------------------------------------
 * Function:   Thread_sum
 * Purpose:    Estimate area using multi-thread
 * In arg:     rank
 * Return val: null
 */
void* Thread_sum(void* rank)
{
    long my_rank = (long) rank;
    long double x;
    long double my_sum = 0.0;
    long double start = a+trap_count*my_rank*h;
    long double end = a+trap_count*(my_rank+1)*h;

    my_sum = (f(start) + f(end))/2.0;

    for (long long i = 1; i <= trap_count; i++)
    {
        while(flag!=my_rank);
        x = start+i*h;
        my_sum += f(x);
        flag = (flag+1)%thread_count;
    }
    sum += my_sum;
    
    return NULL;
}  /* Thread_sum */

我们可以减少临界区代码的执行次数，得到以下的线程函数：

/*------------------------------------------------------------------
 * Function:   Thread_sum
 * Purpose:    Estimate area using multi-thread
 * In arg:     rank
 * Return val: null
 */
void* Thread_sum(void* rank)
{
    long my_rank = (long) rank;
    long double x;
    long double my_sum = 0.0;
    long double start = a+trap_count*my_rank*h;
    long double end = a+trap_count*(my_rank+1)*h;

    my_sum = (f(start) + f(end))/2.0;
    
    for (long long i = 1; i <= trap_count; i++) {
       x = start+i*h;
       my_sum += f(x);
    }

    while(flag!=my_rank);
    sum += my_sum;
    flag = (flag+1)%thread_count;
    
    return NULL;
}  /* Thread_sum */

2. mutexes
   在主函数中，我们对mutexes进行初始化，并且在程序结束前销毁mutexes。
   在之前的主函数中做一些修改：
   增加pthread_mutex_init函数和pthread_mutex_destroy函数：

使用mutexes的线程函数实现如下所示：

/*------------------------------------------------------------------
 * Function:   Thread_sum
 * Purpose:    Estimate area using multi-thread
 * In arg:     rank
 * Return val: null
 */
void* Thread_sum(void* rank)
{
    long my_rank = (long) rank;
    long double x;
    long double my_sum = 0.0;
    long double start = a+trap_count*my_rank*h;
    long double end = a+trap_count*(my_rank+1)*h;

    my_sum = (f(start) + f(end))/2.0;
    for (long long i = 1; i <= trap_count; i++)
    {
        x = start+i*h;
        my_sum += f(x);
    }

    pthread_mutex_lock(&mutex);
    sum += my_sum;
    pthread_mutex_unlock(&mutex);

    return NULL;
}  /* Thread_sum */

3. semaphores
   在主函数中，我们对semaphores进行初始化，并且在程序结束前销毁semaphores。
   在之前的主函数中增加了sem_init和sem_destroy函数：

使用semaphores的线程函数的实现如下所示：

/*------------------------------------------------------------------
 * Function:   Thread_sum
 * Purpose:    Estimate area using multi-thread
 * In arg:     rank
 * Return val: null
 */
void* Thread_sum(void* rank)
{
    long my_rank = (long) rank;
    long double x;
    long double my_sum = 0.0;
    long double start = a+trap_count*my_rank*h;
    long double end = a+trap_count*(my_rank+1)*h;

    my_sum = (f(start) + f(end))/2.0;
    for (long long i = 1; i <= trap_count; i++)
    {
        x = start+i*h;
        my_sum += f(x);
    }

    sem_post(&semaphore_p);
    sum += my_sum;
    sem_wait(&semaphore_p);

    return NULL;
}  /* Thread_sum */

结果

由于结果测试结果众多，所以选取部分结果截图展示，详细结果可见附件的result.txt文件。
当选取4个线程，划分40000个小梯形，计算在f(x)=x²上从0到2的面积。分别使用busy-waiting，mutexes和semaphores三种方式。结果依次如下所示：


由上图我们可以得知，当n比较小的时候，大概小于10的5次方数量级的时候，并行计算的所耗时间和串行计算所耗时间相差不多，准确度也不是很高。当n大于10的5次方数量级时，准确度就越高，并行计算的所耗时间就越短。
而且在相同的线程数量和n的条件下，busy-waiting所花费的时间小于mutexes所花费的时间，小于semaphores所花费的时间。所以说使用semaphores来保护临界区是一个比较理想的做法。

程序运行方法

源代码中包括了pth_trap_busy.c；pth_trap_mutex.c；pth_trap_semap.c；result.txt；timer.h这5个文件，其中result.txt中是运行的控制台结果，timer.h是这3个C文件的头文件，这些都要放在同一个目录下面。
程序运行的环境中是Ubuntu16.04，在控制台下分别输入命令：

gcc -g -Wall -o pth_trap_busy pth_trap_busy.c -lm -lpthread
gcc -g -Wall -o pth_trap_mutex pth_trap_mutex.c -lm -lpthread
gcc -g -Wall -o pth_trap_semap pth_trap_semap.c -lm -lpthread

即可编译以上三个C文件
然后输入以下命令即可运行：

./pth_trap_busy <number of threads> <n> a b
./pth_trap_mutex <number of threads> <n> a b
./pth_trap_semap <number of threads> <n> a b

number of threads是你想运行的线程数；n是所求面积划分出的梯形个数，并且n应该被线程数整除；a应该小于等于比b。
源代码中的注释十分详细，如仍有不懂，可以私信咨询。