共享数据的并发修改

　　因为SRIM主要模拟跟踪一大批入射粒子的运动，每个粒子在靶物质中的运动是相互独立的，所以把每个粒子在靶材料中的作用作为线程并行化，在统计能损过程会遇到如此问题：

　　MTOT[i]=MTOT[i]+E_loss,若把靶在深度方向上的长度均匀划分成N段，则i为第i段靶材料，MTOT[i]用来存储所有粒子在第i段上的能损总量，E_loss是每个粒子在第i段上的能损。此类问题即为浮点运算共享数据的并发修改问题，是并行开发中最常见的问题。对于整数型变量的并发修改我们可利用CUDA提供的原子加法操作，但由于浮点数加法运算的中间结果可能被截断，并非是可结合的，所以在这里不能应用。我们也可以自己书写高级原子操作通过锁定数据结构来实现，但实践证明其在线程数比较多的情况下很耗时。还有一种可行的方法就是，对变量MTOT分配Yn*blockdim*threaddim大小,其中blockdim=(ns_p+threaddim-1)/threaddim即把每个粒子的对应在靶材料中的能损分别存储，然后再并行相加，这种方法在粒子数较小的情况下是可行的，但等到粒子数很大时改变量会占用很大内存。本文提出一种折中的方法，同样给MTOT分配Yn*blockdim*threaddim大小，其中blockdim，threaddim给于较小的定值，如二者都为16。

　　每个线程的起始索引按照如下公式计算int tid=threadIdx.x+blockDim.x*blockIdx.x;

　　当粒子数大于线程数时，在每个线程计算完当前索引任务后接着需要对索引进行递增，其中递增的步长为线程格中正在运行的线程数量。即i=i+blockDim.x*gridDim.x；

　　那么代码框架可写为如下：

　　…
　　#define yn (100)
　　…
　　__global__ void mcc(参数){
　　int tid=threadIdx.x+blockDim.x*blockIdx.x;
　　int i=tid;
　　…
　　while(i<particleNum&&tid<blockDim.x*gridDim.x){
　　…
　　temp=(i%blockDim.x*gridDim.x)*yn+m-1;//m>=0&&m<yn
　　MTOT[temp]=MTOT[temp]+E_loss;
　　…
　　i+= blockDim.x*gridDim.x;
　　}
　　}
　　int main(){
　　int blockdim=16;
　　int threaddim=16;
　　int n_HN_total= yn* blockdim* threaddim;
　　…
　　mcc<<<blockdim,threaddim>>>( d_MTOT，…参数);
　　CUDA_SAFE_CALL(cudaThreadSynchronize());
　　CUDA_ SAFE_CAL L (cudaMemcp y( h_MTOT,d_MTOT,
　　n_HN_total*sizeof(float),cudaMemcpyDeviceToHost));
　　for(int u=0;u<thread_total;u++){
　　for(int v=0;v<yn;v++){
　　MTOT[v]=MTOT[v]+h_MTOT[u*yn+v];
　　…
　　}
　　}
　　…
　　}

　　这种实现方法的好处是：1.无论线程数增加到多大，该代码都是可行的；2.其占用的内存是固定的不会随着线程数的增加而增大。