c语言蒙特卡洛求pi,蒙特卡洛算法求 pi项目(3)

六、并行算法改进

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

6.1减少线程转态转换次数

此方法具体为：在并行程序中使用互斥锁，当某一线程进入临界区后，一次性产生m个随机点，然后再退出临界区，开始对m个点进行计算；与此同时，若另一线程也要进入临界区，则被挂起，等待该线程退出。如此循环，直至两个线程均计算完所要求的点的个数，则计算输出值，程序结束。

算法：

1、确定产生点n的个数和缓冲区m(m<=n)的值，声明互斥锁

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 *)

voidthread(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)<=1)

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)<=1)

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<<31;

//生成伪随机数算法

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<<31;

//生成伪随机数算法

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)<=1)

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)<=1)

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);

}