CUDA C 并行编程
CPU矢量求和运算
/*
//代码4.2.1.1基于CPU的矢量求和运算
//时间：2019.07.20
*/
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#define N 10
void 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];
//在CPU上为数组'a'和'b'赋值
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]);
}
system("pause");
return 0;
}
CPU代码中add()函数通过while循环对矢量a，b中N个元素串行依次加和并存储到矢量c中相应位置。
GPU矢量求和运算
/*
//代码4.2.1.2基于GPU的矢量求和运算
//时间：2019.07.20
*/
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#define N 10
__global__ void add(int *a, int *b, int *c)
{
int tid = blockIdx.x;//计算位于这个索引处的数据
if (tid < N)
c[tid] = a[tid] + b[tid];
}
int main()
{
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));
//在CPU上为数组'a'和‘b’赋值
for (int i = 0; i < N; i++)
{
a[i] = -i;
b[i] = i*i;
}
//将数组'a'和'b'复制到GPU
cudaMemcpy(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复制到CPU
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);
system("pause");
return 0;
}
需要注意的是，上面的代码中使用了一些通用模式：
(1)调用cudaMalloc()在设备上未三个数组分配内存：在其中两个数组(dev_a和dev_b)中包含了输入值，而在数组dev_c中包含了计算结果。
(2)为了避免内存泄漏，在使用完GPU内存后通过cudaFree()释放他们。
(3)通过cudaMemcpy()将输入数据复制到设备中，同时指定参数cudaMemcpyHostToDevice，在计算完成后，将计算结果通过参数cudaMemcpyDeviceToHost复制回主机。
(4)通过尖括号语法，在主机代码main()中执行add()中的设备代码。
尖括号语法N,1说明
<<<N,1>>>
尖括号语法用于指定启动核函数的方式，其中第一个参数表示设备在执行核函数时使用的并行线程块block的数量，第二个参数表示一个block中线程thread的数量。
<<<N,1>>>说明在执行核函数的时候将创建核函数的N个副本，即N个block，每个block中只包含一个thread，并以并行的方式来运行它们。
通用尖括号语法gridsize,blocksize说明
<<<gridsize,blocksize>>>
一个GPU只能运行一个kernal，一个kernal就是一个grid，一个grid中可以包含多个block，一个block可以包含多个thread。
dim3 gridsize(1,2,3)
dim3 blocksize(4,5,6)
<<<gridsize,blocksize>>>
dim3是CUDA提供的用于指明三维结构的结构体类型,上述代码首先指明了grid以三维的形式组织，每一维包含的block数量分别为(1,2,3),接着指明了block也以三维的形式组织，每一维包含的thread数量分别为(4,5,6).
除此之外，dim3还可以以二维的形式组织甚至是一维的形式组织。
dim3 gridsize(7,8)
dim3 blocksize(9)
后面省略掉的维数将自动由编译器用1进行填充。
核函数说明
__global__ void add(int *a, int *b, int *c)
{
int tid = blockIdx.x;//计算位于这个索引处的数据
if (tid < N)
c[tid] = a[tid] + b[tid];
}
其中，最重要的是内置变量blockIdx的使用。
blockIdx内包含三个成员变量，分别是x，y，z，用于指明当前核函数运行在哪个block中，由于上述<<<N,1>>>blocksize只有一维有效，因此单blockIdx.x就足以指明block的索引值。
基于CPU的Julia集
//代码4.2.2.1基于CPU的Julia集
//时间：2019.07.24
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include "cpu_bitmap.h"
#define DIM 1000
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);
}
};
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;
}
void kernal(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 main()
{
CPUBitmap bitmap(DIM, DIM);
unsigned char *ptr = bitmap.get_ptr();
kernal(ptr);
bitmap.display_and_exit();
}
可以看到，Julia图案生成是逐像素的，生成速度很慢。
代码解释：
1.cuComplex
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);
}
};
定义一个结构体用于保存复数值，其中r是实部，i是虚部，magnitude2()计算该复数的模值， operator*()对复数乘法进行重载， operator+()对复数加法进行重载。
2.julia()
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;
}
给定想要计算的图像中的某一个点的位置(x,y),计算该点是否属于julia集合，如果属于则返回1，不属于则返回0.
函数首先将像素坐标转换为复数空间的坐标。为了将复平面的原点定位到图像中心，我们将像素位置移动了DIM/2.为了确保图像的范围为-1.0到1.0,我们将图像的坐标缩放了DIM/2倍。
然后，我们引入了一个scale因数来实现图形的缩放，当前，其被硬编码为1.5。
在计算出复数空间点之后通过计算迭代等式Zn+1=Zn2+C,其中C是任意一复数常量，选择值为-0.8+0.156i恰能够生成一张有趣的图片。
迭代计算200次，每次都判断模值是否大于阈值，如果大于阈值，认为等式是发散的，该点不属于Julia集，返回0，对图片中的这个点不着色。如果200次迭代之后模值仍然小于阈值，认为等式是收敛的，该点属于Julia集合，返回1，对图片中的这个点着色。
3.kernal()
void kernal(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;
}
}
}
核函数对将要绘制的所有点进行迭代，并在每次迭代时调用julia来判断该点是否属于Julia集。如果该点位于集合中，那么函数julia()将返回1，否则将返回0。如果julia()返回1，那么就将点的颜色设置为红色，如果返回0则设置为黑色。
4.main()
int main()
{
CPUBitmap bitmap(DIM, DIM);
unsigned char *ptr = bitmap.get_ptr();
kernal(ptr);
bitmap.display_and_exit();
}
通过工具库创建了一个指定大小的位图图像，并将一个指向位图数据的指针传递给了核函数，核函数填充完这个位图图像之后对其进行显示。
基于GPU的Julia集
//代码4.2.2.2基于GPU的Julia集
//时间：2019.07.24
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include "cpu_bitmap.h"
#define DIM 1000
struct cuComplex{
float r;
float i;
__device__ cuComplex(float a, float b) :r(a), i(b){}//注意，这里的__device__是我自己添加的，指明构造函数也是运行在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);
}
};
__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;
}
__global__ void kernal(unsigned char *ptr)
{
//将threadIdx/blockIdx映射到像素位置
int x = blockIdx.x;
int y = blockIdx.y;
int offset = x + y*gridDim.x;//这里需要注意的是blockIdx.x是二级索引，blockIdx.y是一级索引
//计算这个位置上的值
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 main()
{
CPUBitmap bitmap(DIM, DIM);
unsigned char *dev_bitmap;
cudaMalloc((void **)&dev_bitmap, bitmap.image_size());
dim3 grid(DIM, DIM);
kernal << <grid, 1 >> >(dev_bitmap);
cudaMemcpy(bitmap.get_ptr(), dev_bitmap, bitmap.image_size(), cudaMemcpyDeviceToHost);
bitmap.display_and_exit();
cudaFree(dev_bitmap);
}
需要注意的是：
(1)__device__修饰符表明该函数必须由运行在GPU上的代码调用
(2)struct cuComplex中的构造函数前面也需要添加__device__修饰符，这是因为构造函数也是运行在device上的，原书中没有添加，编译会出错。
(3) blockIdx.x是二级索引，blockIdx.y是一级索引
int x = blockIdx.x;
int y = blockIdx.y;
所以int offset = x + y*gridDim.x;
实验发现，对比CPU版本，GPU版本并没有可观的速度提升。
本章小结
1.dim3 gridsize(1,2,3)
2.dim3 blocksize(4,5)
3,kernal<<<gridsize,blocksize>>>()
4. blockIdx.x [0,0] 范围
blockIdx.y [0,1]
blockIdx.z [0,2]
5. threadIdx.x [0,4]
threadIdx.y [0,3]
6. gridDim.x 1 大小
gridDim.y 2
gridDim.z 3
7. blockDim.x 4
blockDim.y 5
8. int offset = blockIdx.x + blockIdx.y*gridDim.x
blockIdx.x是二级索引，blockIdx.y是一级索引
1138
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
