第九章 CUDA原子函数的合理使用

        本章介绍一类实施原子的函数，简称原子函数。这里的原子与物理学中的原子没有直接关系，但是在概念上有相同之处。从化学的角度上讲，原子是化学物质的基本组成部分，被认为是不可分的。类似的，在CUDA中，一个线程的原子操作可以在不受其他线程的任何操作的影响下完成对某个数据的一套“读-改-写”操作。该套操作也可是说是不可分的。

9.1 完全在GPU中进行归约

        在第八章中几个版本的规约函数中，核函数并没有做完全的计算，而只是将一个长一些的数组d_x变成了一个短一些的数组d_y,后者中每个元素为前者若干元素的和。由两种方法能够在GPU中得到最终结果，一是用另一个核函数将较短的数组进一步归约，得到最终的结果(一个数值)；二是在先前的核函数的末尾利用原子函数进行规约，直接得到最终的结果，本章仅讨论第二种方法。
        本章只介绍第二种方法。为此先考察前一章使用GPU进行规约核函数的最后几行。
        if(tid==0)
{
        d_y[bid]=s_y[0];
}
这个if语句块的作用是将每一个线程块中的归约结果从共享内存s_y[0]复制到全局内存d_y[bid]。为了将不同线程块的部分和s_y[0]累加起来，存放到一个全局内存地址，可以将上述代码改为
if(tid==0)
{
d_y[0]+=s_y[0];
}
遗憾的是，该语句不能实现目的。该语句在每一个线程块的第0号线程都会被执行，但是他们执行的次序是不确定的。在每个线程中，该语句其实可分为两个操作：首先从d_y[0]中读取数据并与s_y[0]相加；然后将结果写入d_y[0],不管次序如何，只有一个线程的“读-写”操作不被其他线程干扰时，才能得到正确的结果，日过一个线程还未将结果写入d_y[0]中，另一个线程就读取了d_y[0],那么这两个线程读取的d_y[0]就一样，这必然会导致错误的结果，考虑bid为0和1的两个线程块中的0号线程，他们将执行如下计算：d_y[0]+=s_y[0];
假如bid为0的线程先读取d_y[0]的值，然后计算d_y[0]+s_y[0],得到了一个数，但是当他还没有来得及将结果写入d_y[0]中时，bid为1的线程也读取了d_y[0]的值。要得到所有线程中的s_y[0]的和，必须使用原子函数：
        if(tid==0)

{
atomicAdd(&d_y[0],s_y[0]);
}
原子函数atomicAdd(adress,val);的第一个参数是待累加变量的地址address中的旧值old读出，计算old+val，然后将计算的值存入地址adress。这些操作在一次原子事务中完成，不会被别的线程中的原子操作干扰。原子函数不能保证各个线程的执行具有特定的顺序，但是能保证每个线程的操作一气呵成，不会被其他线程干扰。这里的if语句能保证每个线程块的s_y[0]只累加一次，去掉if语句会导致结果扩大128倍，且原子操作atomicAdd()的第一个参数是待累加变量的指针，也可以将&d_y[0]写成d_y.
        最后要注意的是，在调用该版本的函数之前，必须将d_y[0]的值初始化为0，在本节中将使用包装函数，对核函数进行调用并做对应的前处理和后处理工作：
 

#include "error.cuh"
#include <stdio.h>

#ifdef USE_DP
    typedef double real;
#else
    typedef float real;
#endif

const int NUM_REPEATS = 100;
const int N = 100000000;
const int M = sizeof(real) * N;
const int BLOCK_SIZE = 128;

void timing(const real *d_x);

int main(void)
{
    real *h_x = (real *) malloc(M);
    for (int n = 0; n < N; ++n)
    {
        h_x[n] = 1.23;
    }
    real *d_x;
    CHECK(cudaMalloc(&d_x, M));
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));

    printf("\nusing atomicAdd:\n");
    timing(d_x);

    free(h_x);
    CHECK(cudaFree(d_x));
    return 0;
}

void __global__ reduce(const real *d_x, real *d_y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int n = bid * blockDim.x + tid;
    extern __shared__ real s_y[];
    s_y[tid] = (n < N) ? d_x[n] : 0.0;
    __syncthreads();

    for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
    {
        if (tid < offset)
        {
            s_y[tid] += s_y[tid + offset];
        }
        __syncthreads();
    }

    if (tid == 0)
    {
        atomicAdd(d_y, s_y[0]);
    }
}

real reduce(const real *d_x)
{
    const int grid_size = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
    const int smem = sizeof(real) * BLOCK_SIZE;

    real h_y[1] = {0};
    real *d_y;
    CHECK(cudaMalloc(&d_y, sizeof(real)));
    CHECK(cudaMemcpy(d_y, h_y, sizeof(real), cudaMemcpyHostToDevice));

    reduce<<<grid_size, BLOCK_SIZE, smem>>>(d_x, d_y, N);

    CHECK(cudaMemcpy(h_y, d_y, sizeof(real), cudaMemcpyDeviceToHost));
    CHECK(cudaFree(d_y));

    return h_y[0];
}

void timing(const real *d_x)
{
    real sum = 0;

    for (int repeat = 0; repeat < NUM_REPEATS; ++repeat)
    {
        cudaEvent_t start, stop;
        CHECK(cudaEventCreate(&start));
        CHECK(cudaEventCreate(&stop));
        CHECK(cudaEventRecord(start));
        cudaEventQuery(start);

        sum = reduce(d_x); 

        CHECK(cudaEventRecord(stop));
        CHECK(cudaEventSynchronize(stop));
        float elapsed_time;
        CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
        printf("Time = %g ms.\n", elapsed_time);

        CHECK(cudaEventDestroy(start));
        CHECK(cudaEventDestroy(stop));
    }

    printf("sum = %f.\n", sum);
}

这里有几点需要值得注意：
        (1)这里的主机数组h_y使用的是栈(stack)内存，而不是堆内存，也就是说传给cudaMemcpy()函数的主机内存可以使栈内存，传输少量数据时安全，但是在传输大量数据时不安全，因为栈的大小是有限的。
        (2)包装函数与核函数同名，这里使用了C++的重载特性：可以同名，但是参数列表可以不一致。原子函数atomicAdd()其实是有返回值的，只不过在上述程序中没有使用。

9.2 原子函数

        原子函数对它的第一个参数指向的数据进行一次读-改-写的原子操作，即一气呵成。第一个参数可以指向全局内存，也可以指向共享内存。对所有参与的线程来说，该“读-改-写”的原子操作是一个线程一个线程轮流做的，但没有固定的次序，另外原子函数没有同步功能。
        从帕斯卡架构开始，在原来的原子函数的基础上引入了两类新的原子函数，例如对于原子函数atomicAdd()来说，从帕斯卡架构起引入了两类新的原子函数，分别是atomicAdd_system()和atomicAdd_block(),前者将原子函数的作用范围扩展到整个同节点的异构系统(包括主机和所有设备)，后者将原子函数的作用范围缩小到一个线程块。
        下面列出所有原子函数的原型，并介绍他么扥功能。我们约定，对每一个线程来说，address所指变量的值在实施与该线程对应的原子函数前为old，在实施对应的原子函数后为new。对每一个原子函数来说，返回值都是old，另外这里介绍的原子函数都是__device__ 函数，只能在核函数中使用。
 

(1)加法：T atomicAdd(T*address,T val);
功能：new=old+val
(2)减法：T atomicSub(T*address,T val);
功能：new=old-val
(3)交换：T atomicExch(T*adress,Tval);
功能：new=val
(4)最小值：T atomicMin(T*address,T val)
new=(old>val)?:val:old
(5)最大值 T atomicMax(T*address,T val)
new=(old>val)?old:val;
(6)T atomicInc(T* address,T val)
功能：new=(old>val)?0:(old+1);
(7)自减 T atomicDec(T*address,T val)
功能：new=((old==0)||(old>val))?val:(old-1)
(8)比较交换(compare and swap)T atomicCAS(T*addres,Tcompare,T val)
功能：new=(old==compare)?val:old.
(9)换位与: T atomicAnd(T *address,T val)
功能：new=old&val
(10)按位或 T atomicOr(T *address,T val);
(11)按位异或 T atomicXor(T*address,T val);
功能：new=old^val

在上面所列的函数中，我们用T表示相关变量的数据类型。在所有原子函数中，atomicCAS()函数是比较特殊的：所有其他原子函数都可以利用它实现。例如，在帕斯卡架构出现以前，atomicAdd()函数不支持双精度浮点数，就有人用atomicCAS()函数自己实现了一个双精度浮点数的atomicAdd(),实现代码如下所示：
 

__device__ double atomicAdd(double* address,double val)
{
unsigned long long *address_as_ull=(unsigned long long *)address;
unsigned long long old=*address_as_ull,assumed;
do
{
    assumed=old;
    old=atomicCAS(address_as_ul,__double_as_longlong(val+__longlong_as_double(assumed)));
}while(assumed!=old);
return __longlong_as_double(old);
}

9.3 例子：邻居列表的建立

        在系统介绍了CUDA中的原子函数之后，用一个建立另据列表的例子进一步说明原子函数的使用。许多计算机模拟都用到了邻居列表，如分子动力学模拟和蒙特卡洛模拟。
        给定N个粒子的坐标ri,如何确定每个粒子的邻居个数及每个邻居的粒子指标呢。其中，两个粒子互相为邻居的判据如下：他们的距离不大于一个给定的横截距离rc。为明确起见，我们以下图所示的原子系统为例进行讨论：

 该图展示的是一个单层多晶石墨烯片中的部分原子坐标。该图中每一个点代表一个碳原子。横坐标x和纵坐标y的单位都是埃，即1e-10m.该体系是二维的，可以明显的看到该图没有显示碳原子之间的化学键。要画出化学键，需要知道哪些原子之间有化学键。这就需要计算邻居列表，即每个原子所有的邻居原子，一个原子和他的每一个邻居原子都有一个碳碳化学键。目的是确定这个体系的邻居列表，从而画出化学键，得到更美观的图。
        我们用一个简单地例子来说明邻居列表的概念。假设有4个粒子，从0-3进行标记。假设他们之间两两之间是邻居，那么每个粒子就有三个邻居，其中第0个粒子的邻居是1，2，3.
        我们的基本算法如下：对于一个给定的粒子，通过比较它与其他粒子的距离来判断相应的粒子对是否为邻居。这是一个计算量正比于N2的算法，对于大粒子群，有一个高效但编程实现相对复杂的O(N)算法。

9.3.1 C++版本的开发

        C++版本的完整程序见以下代码，该代码包含三个部分，1 从文件xy.txt中读取坐标数据；2 调用函数find_neignbor()构建邻居列表并对该函数进行计时。3 将计算出来的邻居列表输出到文件neighbor.txt:
 

#include"error.cuh"
#include<cmath>
#include<iostream>
#include<fstream>
#include<sstream>
#include<string>
#include<vector>

#ifdef USE_DP
typedef double real;
#else
typedef float real;
#endif

int N;//number of atoms
const int NUM_REPEATS=20;//number of timings
const int MN;//maximun number of neibors for each atom
const real cutoff=1.9;//in units of angstrom
const real cutoff_squere=cutoff*cutoff;

void read_xy(std::vector<real>&x,std::vector<real>&y);
void timing(int *MN,int *NL,std::vector<real>x,std::vector<real>y);
void print_neignbor(const int*NN,const int NL);

int main()
{
std::vector<real>x,y;
read_xy(x,y);
N=x.size();
int *NN=(int*)malloc(N*sizeof(int));
int *NL=(int*)malloc(N*MN*sizeof(int));

timing(NN,NL,x,y);
print(NN,NL,x,y);
free(NN);
free(NL);
return 0;
}
void read_xy(std::vector<real>&v_x,std::vector<real>&v_y)
{
    std::ifstream infile("xy.txt");
    std::string line,word;
    if(!inline)
    {
    std::cout<<"Cannot open xy.txt"<<std::endl;
    exit(1);
    }
    while(std::getline(infile,line))
    {
    std::isstringstream words(line);
    if(line.length()==0)
        {
    continue;
        }
    for(int i=0;i<2;i++)
        {
            if(words>>word)
                {
                    if(i==0)
                        {
                            v_x.push_back(std::stod(word));
                        }
                    if(i==1)
                        {
                            v_y.push_back(std::stod(word));
                        }
                }
            else
                {
                    std::cout<<"Error for reading xy.txt"<<std::endl;
                    exit(1);
                }
        }
    }
    infile.close();
}
void find_neighbor(int *NN,int * NL,const real *x,const real *y)
{
    for(int n=0;n<N;n++)
    {
        NN[n]=0;
    }
    for(int n1=0;n1<N;++n1)
    {
        real x1=x[n1];
        real y1=y[n1];
        for(int n2=n1+1;n2<N;++n2)
            {
        real x12=x[n2]-x1;
        real y12=y[n2]-y1;
        real distance_squere=x12*x12+y12*y12;
        if(distance_square<cutoff_square)
                {
                    NL[n1*MN+NN[n1]++]=n2;
                    NL[n2*MN+NN[n2]++]=n1;
                }
            }
    }
}
void timing(int *NN,int *NL,std::vector<real>x,std::vector<real>y)
{
    for(int repeat=0;repeat<NUM_REPEATS;++repeat)
        {
            cudaEvent_t start,stop;
            CHECK(cudaEventCreate(&start));
            CHECK(cudaEventCreate(&stop));
            CHECK(cudaEventRecord(start));
            while(cudaEventQuery(start)!=cudaSuccess){}
            find_neighbor(NN,NL,x.data().y.data());
            
            CHECK(cudaEventRecord(stop));
            CHECK(cudaEventSynchronize(stop));
            float elapsed_time;
            CHECK(cudaEventElapsedTime(&elapsed_time,start,stop));
            std::cout<<"Time="<<elapsed_time<<"ms."<<std<<endl;
            
            CHECK(cudaEventDestroy(start));
            CHECK(cudaEventDestroy(stop));
        }
}

void print_neighbor(const int *NN,const int *NL)
{
    std::ofstream outline("neighbor.txt");
    if(!outfile)
    {
        std::cout<<"Cannot open neighbor.txt"<<std::endl;
        
    }
    for(int n=0;n<N;++n)
    {
        if(NN[n]>MN)
        {
            std::cout<<"Error:MN is too small."<<std::endl;
            exit(1);

        }
        outfile<<NN[n];
        for(int k=0;k<MN;++k)
        {
            if(k<NN[n])
            {
                outfile<<" "<<NL[n*MN+k];
            }
            else
            {
                outfile<<"NaN";
            }
        }
    outfile<<std::endl;
    }
outfile.close();
}

该函数涉及的参数与数据如下：
        (1)总的原子数N.
        (2)每个原子的最多邻居原子数MN，我们这里将其设置为10，代表我们认为不可能有原子拥有十个以上的邻居。实际上，在我们选取的横断面积内，不可能有原子拥有三个以上的邻居。所以将该变量初始化为3，这里需要注意的是，如果该参数取得小了，就有可能出现内存错误；
        (3)判断两个原子是否为邻居的截断距离平方cutoff_square,这里取为1.9E2埃尔的平方。
        (4)整形数组NN。该数组的长度为N，其中NN[n]是第n个粒子的邻居个数。
        (5)整形数组NL。该数组的长度为N*MN，其中NL[n*MN+k]是第n个粒子的第k个邻居的指标。
        (6)单精度或双精度浮点型数组x,y.他们的长度都是N，分别记录每个原子的x坐标和y坐标。

下面对find_neighbor()的实现做一些解释：
        (1)函数的第3-6行将数组NN的所有元素初始化为零，因为后面要对数组元素进行累加。
        (2)从第8行开始，对体系的所有原子进行二重循环。循环变量n1和n2就代表两个可能互为邻居的原子。入股n2是n1的邻居，那么反过来，n1就一定是n2的邻居，从而省去一半的计算。
        (3)在循环体中，首先计算两个原子之间的距离平方distance_square,并与截断距离的平方进行比较，若原子之间的距离平方小于截断距离平方，则给原子n1添加一个邻居原子n2，并将n1的邻居原子数NN[n1]增1，接着给原子n2添加一个邻居原子n1，并将n2的邻居原子数NN[n2]增1，第19行必须使用后置++，这样能保证先利用它的值，而后进行增加。
       用以下命令编译单精度版本
        nvcc -O3 -arch=sm_75 neighbor1cpu.cu
        用以下命令编译双精度版本
        nvcc -O3 -arch=sm_75 -DUSE_DP neighbor1cpu.cu
得到输出文件之后进行绘图得到如下输出结果

 

9.3.2 利用原子操作的CUDA版本

        首先从上面的C++版本的find_neighbor()函数出发，写出一个对应的核函数，以下代码给出了核函数版本find_neighbor_atomic():
 

void __global__ find_neighbor_atomic(int *d_NN,int*d_NL,const real *d_x,const real *d_y,
const int N,const real cutoff_square)
{
int n1=blockIdx.x*blockDim.x+threadIdx.x;
if(n1<N)
{
d_NN[n1]=0;
real x1=d_x[n1];
real y1=d_y[n1];
for(int n2=n1+1;n2<N;++n2)
{
    real x12=d_x[n2]-x1;
    real y12=d_y[n2]-y1;
    real distance_square=x12*x12+y12*y12;
    if(distance_square<cutoff_square&&(n2!=n1))
        {
            d_NL[n1*MN +atomicAdd(&d_NN[n1],1)]=n2;
            d_NL[n2*MN+atomicAdd(&d_NN[n2],1)]=n1;
        }
}
}
}

在使用原子函数时还有一个容易犯的错误，就是没有合理利用返回值。在C++版本的程序中，我们可以将代码
        NL[n1*MN+NN[n1]++]=n2;
        NL[n2*MN*NN[n2]++]=n2;
改写成以下形式
        NL[n1*MN+NN[n1]]=n2;
        NL[n1]++;
        NL[n2*MN+NN[n2]]=n1;
        NN[n2]++;
他们都是正确的。然而如果类似地使用原子函数的代码
d_NL[n1*MN+d_NN[n1]]=n2;
atomicAdd(&d_NN[n1],1);
d_NL[n2*MN+d_NN[n2]]=n1;
atomicAdd(&d_NN[n2],1);
那就大错特错了。这是因为，改写后第一行使用的d_NN[n1]的值并不能保证的第二行的原子函数读取的d_NN[n1]的值，而有可能是被别的线程的原子函数改动过得值。这是在使用原子函数所访问变量旧值时需要注意的一点，下面是一种正确的写法：
        int tmp1=atomicAdd(&d_NN[n1],1);
d_NL[n1*MN+tmp1]=n2;
int tmp2=atomicAdd(&d_NN[n2],1);
d_NL[n2*MN+tmp2]=n1;
在这种写法中，用临时变量tmp1和tmp2保留了d_NN[n1]和d_NN[n2]的旧值。

9.3.3 不用原子操作的CUDA版本

        对于一个使用了原子函数的问题，有时也可以通过改变算法，去掉对原子函数的使用。在上面的问题中，之所以需要使用原子函数，是因为我们省去了一半对距离的判断，使得不同线程间可能同时写入同一个全局内存地址。如果我们不省去那一半对距离的判断，就可以在一个线程中仅仅写入与之对应的全局内存地址，从而不需要使用原子函数，所写出来的代码如下所示：
 

void __global__ find_neighbor_no_atomic(
int *d_NN,int *d_NL,const real *d_x,const real *d_y,
const int N,const rel cutoff_squere
)
{
    int n1=blockIdx.x*blockDim.x+threadIdx.x;
    if(n1<N)
        {
            int count=0;
            real x1=d_N[n1];
            real y1=d_y[n1];
           
            for(int n2=0;n2<N;++n2)
                {
real x12=d_x[n2]-x1;
real y12=d_y[n2]-y1;
 real distance_squere=x12*x12+y12*y12;
if((distance_squere<cutoff_square)&&(n2!=n1))
                    d_NL[(count++)*N+n1]=n2;
                }
        }
        d_NN[n1]=count;
}

相比于使用原子函数的版本，主要由如下几点变化：
        (1)第13行中的循环变量n2的值从0开始，而不是从n1+1开始。也就是说，对每个粒子n1，我们逐一考察所有其他粒子是不是他的邻居。
        (2)正因为循环变量n2的值包含n1，所以我们在第18行需要判断n2是否等于n1.这里，我们使用了一个合并判断语句的技巧，因为n2!=n1且distance_squere<cutoff_square同时成立时，才认为n2时n1的邻居。
        在该复合判断中，还有一个值得注意的地方是，我们将距离的判断放在&&的前面，将粒子身份放在&&的前面，这是因为对于距离的判断得到假的概率要小得多，而逻辑与操作在得知前面的判断结果为假之后就不会继续判断了，可自行修改探测时间。
        (3)另外一个值得一提的优化策略是用一个寄存器变量count减少了内存变量d_NN[n1]的访问。
        (4)最后改变了邻居列表的排列方式：将d-NL[n1*MN+count++]改成了d-NL[(count++)*N+n1]这是因为n1的变化步调与threadIdx.x一致，这样修改之后，对全局内存d_NL的访问将是合并的。