GPU高性能编程CUDA实战-第四章

//基于GPU的矢量求和
#define N 10

__global__ void add(int* a, int* b, int* c) {
    int tid = blockIdx.x;    //blockIdx是一个内置变量，在CUDA运行时已经预先定义了这个变量
    if (tid < N)    //检查tid是否小于N，是为了防止内存非法访问。如果发生这类问题，HANDLE_ERROR()宏将检测出这种情况并发出警告
        c[tid] = a[tid] + b[tid];
}

int main(void) {
    int a[N], b[N], c[N];
    int *dev_a, *dev_b, *dev_c;

//在GPU上分配内存
    HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));

//在CPU上为数组'a'和'b'赋值
    for (int i = 0; i < N; i++) {
        a[i] = -i;
        b[i] - i * i;
    }

//将数组'a'和'b复制到GPU'
    HANDLE_ERROR(cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
    HANDLE_ERROR(cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));

//通过尖括号语法，在主机代码main()中执行add()中的设备代码
//N表示设备在执行核函数时使用的并行线程块数量，并行线程块集合称为一个线程格
    add<<<N, 1 >>> (dev_a, dev_b, dev_c);    

//将数组'c'从GPU复制到CPU
    HANDLE_ERROR(cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost));


//显示结果
    for (int i = 0; i < N; i++) {
        printf("%d+%d=%d\n", a[i], b[i], c[i]);
    }

//释放在GPU上分配的内存
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return 0;
}
//在启动线程块数组时，数组每一维的最大数量都不能超过65535，这是一种硬件限制。

迭代式：

//基于CPU的Julia集
//main函数通过工具库创建了一个大小合适的位图图像，并将一个指向位图数据的指针传递给核函数
int main(void) {    
    CPUBitmap bitmap(DIM, DIM);    //创建一个DIM*DIM大小的位图图像
    unsigned char* ptr = bitmap.get_ptr();

    kernel(ptr);

    bitmap.diaplay_and_exit();
}

//核函数对将要绘制的所有点进行迭代，并在每次迭代时调用julia()来判断该点是否属于julia集
//如果该点位于集合中，函数julia()将返回1，点的颜色设置成红色；否则返回0，并设置成黑色
void kernel(unsigned char* ptr) {
    for (int y = 0; y < DIM; y++) {
        for (int x = 0; x < DIM; x++) {
            int offset = x +y * DIM;

            int juliaValue = julia(x, y);
            ptr[offset * 4 + 0] = 255 * juliaValue;
            ptr[offset * 4 + 1] = 0;
            ptr[offset * 4 + 2] = 0;
            ptr[offset * 4 + 3] = 255;
        }
    }
}

int julia(int x, int y) {
    const float scale = 1.5;    //引入scale因数来实现图形的缩放

    //将像素坐标转换为复数空间坐标
    float jx = scale * (float)(DIM / 2 - x) / (DIM / 2);    //为了将复平面的原点定位到图像中心，将像素位置移动了DIM/2
    float jy = scale * (float)(DIM / 2 - y) / (DIM / 2);    //为了确保图像范围为-1.0到1.0，将图像坐标缩放了DIM/2倍

    //在计算出复空间中的点后，需要判断这个点是否属于Julia集，选择的值是-0.8+0.156i,如果想看其他的Julia集，可以修改这个常量
    cuComplex c(-0.8, 0.156);
    cuComplex a(jx, jy);

    //共有两百次迭代，每次迭代计算完成后，都会判断结果是否超过某个阈值（这里是1000）
    //如果超过了1000，那么等式就是发散的，将返回0；反之，如果所有200次都计算完毕后，结果仍小于1000，就认为这个点属于该集合，并给调用者kernel()返回1
    int i = 0;
    for (i = 0; i < 200; i++) {
        a = a * a + c;
        if (a.magnitude2() > 1000)
            return 0;
    }
    return 1;
}

//定义一个通用结构来保存复数值，包含两个成员：单精度的实部r和单精度的虚部i
struct cuComplex {
    float r;
    float i;
    cuComplex(float,float b):r(a),i(b){ }
    float magnitude2(void) { return r * r + i * i; }
    cuComplex operator*(const cuComplex& a) {
        return cuComplex(r * a.r - i * a.i, i * a.r + r * a.i);
    }
    cuComplex operator+(const cuComplex& a) {
        return cuComplex(r + a.r, i + a.i);
    }
};

//下面这段没跑出来...

//基于GPU的Julia集
#define DIM 1000
int mian(void) {
    CPUBitmap bitmap(DIM, DIM);   //创建一个DIM*DIM大小的位图图像
    unsigned char* dev_bitmap;    //声明一个指针dev_bitmap,用来保存设备上数据的副本

    HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));   //为了保存数据，通过cudaMalloc()来分配内存

    dim3 grid(DIM, DIM);    //dim3表示一个三维数组，可用于指定启动的线程块数量，最后一维大小为1
    kernel << <grid, 1 >> > (dev_bitmap);    //将dim3变量grid传递给CUDA运行

    HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(),    //将结果复制回主机
        dev_bitmap,
        bitmap.image_size(),
        cudaMemcpyDeviceToHost));
    bitmap.display_and_exit();

    cudaFree(dev_bitmap);
}

//CPU与GPU版本之间的差异之一在于kernel()的实现
__global__ void kernel(unsigned char* ptr) {    //将kernel()声明为一个__global__函数
    //将threadIdx/BlockIdx映射到像素位置,在(0,0)和(DIM-1,DIM-1)之间的每个像素点(x,y)都能获得一个线程块
    int x = blockIdx.x;
    int y = blockIdx.y;
    int offset = x + y * gridDim.x;    //输出缓冲区ptr中的线性偏移，通过内置变量gridDim来计算

    //现在计算这个位置上的值
    int juliaValue = julia(x, y);
    ptr[offset * 4 + 0] = 255 * juliaValue;
    ptr[offset * 4 + 1] = 0;
    ptr[offset * 4 + 2] = 0;
    ptr[offset * 4 + 3] = 255;
}

//判断某点是否属于Julia集
__device__ int julia(int x, int y) {
    const float scale = 1.5;    //让图像更大（1.5倍）
    float jx = scale * (float)(DIM / 2 - x) / (DIM / 2);    //让图像开始时位于窗口中心
    float jy = scale * (float)(DIM / 2 - y) / (DIM / 2);

    cuComplex c(-0.8, 0.156);
    cuComplex a(jx, jy);

    //执行迭代运算，判断该像素点是否属于Julia集，是的返回1
    int i = 0;
    for (i = 0; i < 200; i++) {
        a = a * a + c;
        if (a.magnitude2() > 1000)
            return 0;
    }
    return 1;
}

//定义cuComplex结构，用单精度浮点数来保存复数的实部和虚部
struct cuComplex {
    float r;
    float i;
    __device__ cuComplex(float a, float b) :r(a), i(b) { }  //书上代码少了__device__
    __device__ float magnitude2(void) {
        return r * r + i * i;   //求复数的模
    }
    __device__ cuComplex operator*(const cuComplex& a) {
        return cuComplex(r * a.r - i * a.i, i * a.r + r * a.i);
    }
    __device__ cuComplex operator+(const cuComplex& a) {
        return cuComplex(r + a.r, i * a.i);
    }
};