并行Pi的计算

主要考虑两种PI的计算方法 ， 一种为积分法的用共享内存试的GPU实现， 另外一种为在面积估算MPI伪代码实现。

1.积分法求pi

// 计算PI 
// 算法：在[0,1]范围内积分 1/(1+x*x)

// 做积分，把每一个block中的和放入sum对应的下标中
__global__ void kernel1(float* sum, int num){
	int gid = blockIdx.x*blockDim.x+threadIdx.x;//线程索引
	
	float temp;
	extern float __shared__ pi[]; // 每一个block中的共享内存
	pi[threadIdx.x] = 0.0; // 初始化
	__syncthreads();

	// 计算积分 ， 放入对应线程里
	while(gid < num){
		temp = (gid + 0.5f) / num;
		pi[threadIdx.x] += 4.0/(1 + temp*temp);
		gid += blockDim.x * gridDim.x
	}

	// 对每一个block归约求和
	for(int i=(blockDim.x>>1);i>0;i>>=1){
       if(threadIdx.x<i){
          pi[threadIdx.x] += pi[threadIdx.x+i];
       }
       __syncthreads();
    }
    if(threadIdx.x==0)
       sum[blockIdx.x]=s_pi[0];
    }
}

// 对sum在执行一次归约求和

__global__ void  kernel2(float *sum,int num,float *pi){
   int id=threadIdx.x;
   extern float __shared__ s_sum[];
   s_sum[id] = sum[id];
    __syncthreads();
    for(int i=(blockDim.x>>1);i>0;i>>=1){
       if(id<i)
       s_sum[id]+=s_sum[id+i];
       __syncthreads();
}
	if(id==0){
    	*d_pi=s_sum[0]/num;
	}
 }

2.面积法求pi

/*  算法
1.Inscribe a circle in a square
2.Randomly generate points in the square
3.Determine the number of points in the square that are also in the circle
Let r be the number of points in the circle divided by the number of points in the square
4.PI ~ 4*r   
*/
npoints = 10000
circle_count = 0

p = number of tasks
num = npoints/p

find out if I am MASTER or WORKER 

do j = 1,num 
  generate 2 random numbers between 0 and 1
  xcoordinate = random1
  ycoordinate = random2
  if (xcoordinate, ycoordinate) inside circle
  then circle_count = circle_count + 1
end do

if I am MASTER

  receive from WORKERS their circle_counts
  compute PI (use MASTER and WORKER calculations)

else if I am WORKER

  send to MASTER circle_count

endif

两种计算模式， 但其实这两种计算方法都比较适合在GPU上实现

详细的过程可以参考 ： 风晨的CUDA入门 和 Introduction to Parallel Computing