parallel programming in CUDA C（GPU并行程序实现数组求和 Julia set）

前言

我们这节会学习到：
Ⅰ.CUDA在实现并行性时采用的一种重要方式
Ⅱ.用CUDA C编写第一段并行代码

一、Summing vector

#define N 10void add(int *a, int *b, int *c){int tid = 0; //这是第0个CPU，因此索引从0开始while(tid<N){c[tid] = a[tid] + b[tid];tid += 1; // 由于只有一个CPU，因此每次递增1}
}int main(){int a[N],b[N],c[N];for(int i = 0; i < N; ++i){a[i] = -i;b[i] = i * i;}add(a, b. c);for(int i = 0; i < N; ++i){printf("%d + %d = %d\n", a[i], b[i], c[i]);}return 0;
}

这个很简单，就是一个循环求和：把数组a和b保存至c中。
啥意思嘛，不是说并行程序吗？
这是写一个C语言的简单的数组求和，我们以这个为底版，改造成一个GPU的并行版本哦。大家可以想一下怎么写，然后继续往下看哦。

二、GPU vector sums

先看一下，main函数

##define N 10
#include <stdio.h>int main(void) {int a[N], b[N], c[N];int* dev_a, * dev_b, * dev_c;//在GPU上分配内存cudaMalloc((void**)&dev_a, N * sizeof(int));cudaMalloc((void**)&dev_b, N * sizeof(int));cudaMalloc((void**)&dev_c, N * sizeof(int));for (int i = 0; i < N; i++) {a[i] = -i;b[i] = i * i;}//将数组a，b复制到GPUcudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice);cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice);add << <N, 1 >> > (dev_a, dev_b, dev_c);//将数组c从GPU复制到CPUcudaMemcpy(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;
}

add函数

__global__ void add(int *a, int *b, int *c){int tid = blockIdx.x;  //该索引的数据if(tid < N)c[tid] = a[tid] + b[tid];
}

大家是不是注意到了，这次调用add是这样了：
add<<<N,1>>>(dev_a, dev_b, dev_c);
之前是 <<<1,1>>>
现在就为大家解密一个！
这里的第一个参数是代表在GPU的kernel中执行的parallel blocks(并行线程块) 的数目，我们把这些并行线程块叫做grid。
这就带来一个问题，我们如何知道现在运行的线程块是哪个呢？
int tid = blockIdx.x; //该索引的数据
大家看一下这里！
没看到定义啊？这是一个CUDA runtime内置的变量。
那为什么是x呢？大家可以去官网查一下，我们就用到x，其他的用不到哦。
然后GPU中会有N个add副本，tid的值从0 ~ N-1（就是blockIdx.x，也就是这里<<<N,1>>>设定的），分别执行。
if(tid < N) 如果没有这个条件，和C语言会有相似的结果：内存越界。
在这里插入图片描述

三、A fun example（Julia Set）

朱利亚集合可以由下式进行反复迭代得到：
对于固定的复数c，取某一z值（如z = z0），可以得到序列
这一序列可能反散于无穷大或始终处于某一范围之内并收敛于某一值。我们将使其不扩散的z值的集合称为朱利亚集合。
这是百度的，我们要用并行程序实现它。
老规矩先来一个CPU的版本：

1.CPU的Julia集

代码如下（示例）：

#define DIM 1000int main( void ) {CPUBitmap bitmap( DIM, DIM );unsigned char *ptr = bitmap.get_ptr();kernel( ptr );bitmap.display_and_exit();
}

我们创建了一个 bitmap，然后把它的指针传到kernel函数中。

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 );// 像素数据的存储格式是 RGBA,对于每一个像素，其包含了 4 个字节分别对应 RGBA 这 4 个通道的数据,大家可以改变一下试试ptr[offset*4 + 0] = 255 * juliaValue;ptr[offset*4 + 1] = 0;ptr[offset*4 + 2] = 0;ptr[offset*4 + 3] = 255;}}}

遍历所有的点，掉用 julia 进行判断：属于集合返回1（红色），不属于返回0（黑色）。

int julia( int x, int y ) { const float scale = 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);int i = 0;for (i=0; i<200; i++) {a = a * a + c;if (a.magnitude2() > 1000)return 0;}return 1;
}

我们首先将像素坐标转换为复数空间的坐标，为了将复平面的原点定位到图像中心，代码将像素位置移动了DIM/2（我们的图像范围是（DIM，DIM）），然后，为了确保图像的范围为-1.0到1.0，我们将图像的坐标缩放了DIM/2倍。就是说原来的点是(x,y)，转换之后是( (DIM/2-x)/(DIM/2), (DIM/2-y) / (DIM/2) )。在计算处复空间中的点之后，需要判断这个点是否属于Julia集。通过迭代判断（本示例迭代200次，在每次迭代完成后，都会判断结果是否超过某个阈值），如果属于集合，就返回1，否则，返回0。
scale
从生成图像（比如最终想可视化出 Julia 集的样子）的角度来看，不同的scale值会影响图像的缩放程度以及细节呈现情况。较大的scale值可能会让图像看起来被拉伸、放大，能够展现出更多局部的细节特征（但也可能导致图像超出显示范围等问题，需要配合其他参数调整）；较小的scale值则会让图像整体收缩，看到的是更宏观的样子，可能会丢失一些细节信息。
这个 cuComplex c(-0.8, 0.156); 是我们设定的，这样会生成的图片很有趣，大家自己改动一下，看看别的效果怎么样。

struct cuComplex {float   r;float   i;cuComplex( float a, 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);}
};

2.GPU的Julia集

int main( void ) {DataBlock   data;CPUBitmap bitmap( DIM, DIM, &data );unsigned char    *dev_bitmap;cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() );data.dev_bitmap = dev_bitmap;dim3    grid(DIM,DIM);kernel<<<grid,1>>>( dev_bitmap );cudaMemcpy( bitmap.get_ptr(), dev_bitmap,bitmap.image_size(),cudaMemcpyDeviceToHost );cudaFree( dev_bitmap );bitmap.display_and_exit();
}

我们最关心的是有多少并行线程运行在kernel，因为每个点都是独立计算，互不影响，这样就可以开启和我们所计算的点一样多的线程
运行在GPU的函数。我们提到过 blockIdx.x，但这次是坐标系，我们得用一个二维的表示方法：
dim3 grid(DIM,DIM);
类型dim3并不是标准C定义的类型，它可以表是一个三维数组，至于为什么不直接用二维数组，CUDA开发人员主要是为了日后的扩展，所以用三维数组来表示二维数组，数组的第三维默认为1。下面的代码将线程块grid传递给CUDA运行时：
kernel<<<grid,1>>>(dev_bitmap);

__global__ void kernel( unsigned char *ptr ) {// map from blockIdx to pixel positionint x = blockIdx.x;int y = blockIdx.y;int offset = x + y * gridDim.x;// now calculate the value at that positionint 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;
}

__device__ int julia( int x, int y ) {const float scale = 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);int i = 0;for (i=0; i<200; i++) {a = a * a + c;if (a.magnitude2() > 1000)return 0;}return 1;
}

struct cuComplex {float   r;float   i;// 这里要给给它的构造函数前面加个 __device__。问题就是 cuComplex 类的构造函数被定义为了一个 host 函数，然后在 device 函数 julia 上被调用了。本书的代码没有加，也许是以前版本的 CUDA 是允许这种行为？__device__ cuComplex( float a, float b ) : r(a), i(b)  {}__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);}
};