//基于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);
}
};