蒙特卡洛算法求 pi项目(3)

xuyuanchao

六、 并行算法改进
为了改进并行算法，得到更高的加速比，有两种途径可以尝试：减少线程状态转化次数和使用可并行的随机数产生算法。简介如下：

6.1 减少线程转态转换次数
       此方法具体为：在并行程序中使用互斥锁，当某一线程进入临界区后，一次性产生m个随机点，然后再退出临界区，开始对m个点进行计算；与此同时，若另一线程也要进入临界区，则被挂起，等待该线程退出。如此循环，直至两个线程均计算完所要求的点的个数，则计算输出 值，程序结束。
      
算法：
       1、确定产生点n的个数和缓冲区m（m）的值，声明互斥锁
2、某一线程进入临界区，上锁
3、该线程一次性生成m个数，其他线程若想进入则挂起等待
4、该线程退出临界区，解锁 ，开始对刚才生成的随机点进行计算
5、重复2-4步，直至每个线程均完成对所要求点的操作
6、统计COUNTi的值
7、计算 的值

在此算法中，每一线程因为争用rand（）函数而产生的状态转化次数范围为[0， ]，平均次数为 ，调整m的值，使生成m个随机点的时间与对m个随机点进行计算的时间相等时，则算法执行速度可达到最大值，即加速比最大。
示例程序参见附件Pmutex.c。


//并行算法

#include
#include
#include
#include
#include
#include


long cs=0; //循环次数
long count=0; //主线程有效次数
long count_thread=0; //thread线程有效次数
struct timeval start, finish; //定义开始结束时间
double diffsec,diffusec;
long t;//每次生成数据数量
pthread_mutex_t mutex;
long double *data_thread,*data_main;

//void *thread(void *)
void  thread(void)
{
   
    int i=0,j=0;
    double x,y;
    //long double *data_thread;
    long double data_thread[t];
    for(i=0;i
    {
    //  data_thread=new (long double)[t];
        pthread_mutex_lock (&mutex);//lock
        for(j=0;j
            data_thread[j]=(long double)rand()/(long double)RAND_MAX;
        pthread_mutex_unlock(&mutex);//unlock
        for(j=0;j
        {
            x=data_thread[j];
            y=data_thread[j+1];
            if((x*x+y*y)
                count_thread++;
        }
        //delete []data_thread;
    }
}

//主线程计算量为总数的一半
int main(void)
{
    printf("Please input the number:");
    scanf("%d",&cs);
        cs=cs*1000000;
    printf("Please input the number for generating the data once:");
    scanf("%d",&t);
    t=t*100000;
    pthread_t id;
    int ret;
    srand( (unsigned)time( NULL ) );
    pthread_mutex_init(&mutex,NULL);//声明互斥锁
    ret=pthread_create(&id,NULL,(void *)thread,NULL); //创建thread线程
    gettimeofday(&start,NULL); //记录开始时间
    int i=0,j=0;
    double x,y;
    long double data_main[t];   
    for(i=0;i
    {
    //  data_main=new (long double)[t];
        pthread_mutex_lock(&mutex);//lock
        for(j=0;j
            data_main[j]=(long double)rand()/(long double)RAND_MAX;
        pthread_mutex_unlock(&mutex);//unlock
        for(j=0;j
        {
            x=data_main[j];
            y=data_main[j+1];
            if((x*x+y*y)
                count++;
        }
        //delete []data_main;
    }
    pthread_join(id,NULL); //两线程同步
    count+=count_thread;
    gettimeofday(&finish,NULL); //记录结束时间
    //计算所耗时间
    diffsec=(double)finish.tv_sec-start.tv_sec;
    diffusec=(double)(finish.tv_usec-start.tv_usec)/(double)1000000;
    diffsec+=diffusec;
    //输出结果
    printf("Cost time=%f second \n",diffsec);
    printf("roop times=%d \n",cs);
    printf("effective times=%d \n",count);
    printf("pi= %f \n",4*(double)count/(double)cs);
    return(0);
}


6.2 使用可并行的随机数生成函数
       生成随机数最常用的方法为线性同余法，其C语言源代码如下：
       //myrand()用的种子
unsigned static Y =568731;
unsigned d=1
//生成伪随机数算法
double inline myrand()
{   
        Y=(15625*Y+22221)%d;
        return (double)Y/(double)(d-1);
}
通过改变种子的值，算法可生成不同的伪随机数列并且可以满足多个处理器同时调用。但调用所需时间略大于调用系统库函数rand（）。(调用myrand()函数的串行算法，见附件Smyrand.c)
示例程序见附件Pmyrand.c

#include
#include
#include
#include
#include
#include
long cs=0; //总循环次数
long count=0; //主线程有效次数
long count_thread=0; //thread线程有效次数
struct timeval start, finish; //定义开始结束时间
double diffsec,diffusec;
long cs_thread=0,cs_main=0;

//myrand()用的种子
unsigned static Y =568731;
unsigned d=1
//生成伪随机数算法
double inline myrand()
{   
    Y=(15625*Y+22221)%d;
    return (double)Y/(double)(d-1);
}

//thread线程计算量为总数的0.4
void thread(void)
{
   
    int i=0;
    double x,y;
    for(i=0;i
    {
        
        x=myrand();
        y=myrand();
        if((x*x+y*y)
            count_thread++;
    }
}

//主线程计算量为总数的0.6
int main (void)
{
    printf("Please input the number:");
    scanf("%d",&cs);
    cs=cs*1000000;
    cs_main=cs*3/5;
    cs_thread=cs*2/5;
    pthread_t id;
    int ret;
    srand( (unsigned)time( NULL ) );
    ret=pthread_create(&id,NULL,(void *) thread,NULL); //创建thread线程
    gettimeofday(&start,NULL); //记录开始时间   
    int i=0;
    double x,y;
    for(i=0;i
    {      
        x=(long double)rand()/(long double)RAND_MAX;
        y=(long double)rand()/(long double)RAND_MAX;   
        if((x*x+y*y)
            count++;
    }
    pthread_join(id,NULL); //两线程同步
    count+=count_thread;
    gettimeofday(&finish,NULL); //记录结束时间   
    //计算时间差
    diffsec=(double)finish.tv_sec-start.tv_sec;
    diffusec=(double)(finish.tv_usec-start.tv_usec)/(double)1000000;
    diffsec+=diffusec;
    //输出结果
    printf("Cost time=%f second \n",diffsec);  
    printf("roop times=%d \n",cs);
    printf("effective times=%d \n",count);
    printf("pi= %f \n",4*(double)count/(double)cs);
    return(0);
}




本文来自ChinaUnix博客，如果查看原文请点：http://blog.chinaunix.net/u2/86537/showart_1964955.html