CUDA学习与尝试:Hello,World

学习一个东西总要记下来，最近联系了实验室准备做有关等体的计算，浩浩的物理系非资深码农生涯也将进入正轨。万事开头难，当然从构架走起，开始学习CUDA，第一步使用CUDA重构以前做过的建议科学运算算例，诸如有限差分解热流啦，快速傅里叶变换啦，多体运动模拟啦，哈密顿原理模拟啊，以及完成计算物理作业。然后在CFD的实践中全部采用CUDA构架。

目录

CUDA的一些扯淡

新入手的电脑已经达到了睿频4 cores 3.5Ghz的能力，作为八线程出现已经相当的强大，以前也尝试写过一些多线程的代码，代表是

使用Java做刚体力学模拟，然后使用遗传算法用刚体力学模拟结果作为控制器的评价值来演化PIDs的控制器，说白了也就是是一个16线程的大玩具。没有什么实际效果。

不过毕竟是并行的玩意，相比下左图所示的线程池竟然要用效率低下的人脑抢占式c(brain)pu时间分片来完成，而且在线程切换时候效率极其低下，cpu在线程切换时候面临着效率极低的窘境，而且常常因为过热，过冷，能量不足，能量过剩而罢工，甚至会因为视觉系统受到特定人物的刺激而开始大量的宕机（发呆）。不得不喟叹人生之艰难，恩，第一部分因为和颜色相关的理由涉及到个人隐私，作为本期期间最大的支线任务之一就先隐去了。

扯了这么多，几百个线程跑代码确实是一件很爽的事情，是想如果我有分身，一个去做作业一个去复习考试一个去学CFD一个去学CUDA，然后本体去恩（zhui）你（mou）懂（gu）得（niang），那么人生应该如何的美好。不过事情总得一件件来我也没有分身。现在就从CUDA的hello,world学起，讲起。

学习资源：

http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html

我这边并不要做一个翻译工作，而是在学习的过程中讲出来，讲出来自己更懂，也能分享知识。

Hello,World以及线程模型

先看官方给出的一个样例：

1. // Kernel definition
2. __global__ void VecAdd(float* A, float* B, float* C)
3. {
4.     int i = threadIdx.x;
5.     C[i] = A[i] + B[i];
6. }
7.  
8. int main()
9. {
10.     …
11.     // Kernel invocation with N threads
12.     VecAdd<<<1, N>>>(A, B, C);
13.     …
14. }

十分简单，这是一个计算长矢量相加的代码。十分有效的露出了模型的真面目：

根据目前的观察CUDA基于一种非常接近科学运算的模型：块-线程模型，

每个块有若干线程（1024），可以有很多很多块，都可以使用一个三维向量来定位。

而在本例中，长矢量使用threadIdx.x来定位，这也就是cuda模型的漂亮之处，只需要一个三维向量就可以确定每个模块的位置，这个理念虽然和普通意义上的编程相去甚远，但是对于科学运算（大量的矩阵运算）是非常漂亮的，而且在大多数时候矩阵止于二维三维即可，因为用cuda做模拟的化空间也不过三维。

至于具体的线程运行顺序，那是NVIDIA应该考虑的事情，不过根据guide中的介绍，其实和c的多维数组储存没什么区别。

恩，现在我们有了一个模型和一个写好的代码，如何编译它并且运行呢，我们以如下环境为例

Macbook pro retina

System : OS X 10.9

CPU: i7-4850HQ

Memory  16 GB 1600 MHz DDR3

Graphics  NVIDIA GeForce GT 750M 2048 MB

编译环境：

gcc (MacPorts gcc49 4.9-20131027_0) 4.9.0 20131027 (experimental),GNU Make 3.81

开发工具：vim

先自己去装toolkit。参照

http://docs.nvidia.com/cuda/cuda-getting-started-guide-for-mac-os-x/index.html

配置~/.bash_profile，加入

1. export PATH=/Developer/NVIDIA/CUDA–5.5/bin:$PATH
2. export DYLD_LIBRARY_PATH=/Developer/NVIDIA/CUDA–5.5/lib:$DYLD_LIBRARY_PATH

重新打开终端或者重启电脑，然后

1. ~$ nvcc —version
2. nvcc: NVIDIA (R) Cuda compiler driver
3. Copyright (c) 2005–2013 NVIDIA Corporation
4. Built on Thu_Sep__5_10:17:14_PDT_2013

Cuda compilation tools, release 5.5, V5.5.0

恩，安装完成了，于是我们来hello个world

把代码写入(vim 有cuda的高亮，赞一个)

1. //test.cu
2. #include“stdio.h”
3. #include<cuda_runtime.h>
4. #define N 100
5. // Kernel definition
6. __global__ void VecAdd(float* A, float* B, float* C)
7. {
8.     int i = threadIdx.x;
9.     C[i] = A[i] + B[i];
10. }
11.  
12. int main()
13. {
14.     // Kernel invocation with N threads
15.     printf(“Hello,World\n”);
16.     float *A=new float[100],*B=new float[100],*C=new float[100];
17.     for(int i=0;i<N;i++)
18.     {
19.         A[i]=i;
20.         B[i]=i*i;
21.     }
22.     VecAdd<<<1, N>>>(A, B, C);
23.     for(int i=0;i<N;i++)
24.     {
25.         printf(“%f “,C[i]);
26.     }
27.     printf(“\n”);
28. }

以及makefile

1. test:test.o
2.     nvcc –o test test.o
3. test.o:test.cu
4.     nvcc –I /Developer/NVIDIA/CUDA–5.5/include –c test.cu
5. clean:
6.     rm *.o test

make，潇洒通过

运行，

1. CUDA$ ./test
2. Hello,World
3. 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 …………

咦。。完全没有运行。。。似乎是内存管理的问题

注意这样一个问题：刚才的程序完全是把电脑里面的主ram撇给gpu跑，这不是扯淡嘛！

其实cuda的gpu内存分为几种，全局内存，共享内存，线程独占内存，但是无论如何。。。host设备的内存当然是不能跑的,

GPU线性存储器使用 cudaMalloc()分配,通过 cudaFree()释放,使用 cudaMemcpy()在设备和主机间传输

–摘自GUIDE book，风辰翻译

所以，无论如何，先把数据跑到GPU里面吧！

于是修改代码如下

1. __global__ void VecAdd(float* A, float* B, float* C)
2. {
3.     int i = threadIdx.x;
4.     for(int j=0;j<1000;j++)
5.         C[i] = (A[i] * B[i]);
6. }
7.  
8. void cpu_VecAdd(int i,float* A, float* B, float* C)
9. {
10.     for(int j=0;j<1000;j++)
11.         C[i] = (A[i] * B[i]);
12. }
13. int main()
14. {
15.     // Kernel invocation with N threads
16.     printf(“Hello,World\n”);
17.     float *A=new float[N],*B=new float[N],*C=new float[N];
18.     for(int i=0;i<N;i++)
19.     {
20.         A[i]=i;
21.         B[i]=2*i;
22.     }
23.  
24.     size_t size = N * sizeof(float);
25.  
26.     float* d_A;
27.     cudaMalloc(>d_A, size);
28.     float* d_B;
29.     cudaMalloc(>d_B, size);
30.     float* d_C;
31.     cudaMalloc(>d_C, size);
32.     float *e=new float[N];
33.     long st=getCurrentTime();
34.     cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);
35.     cudaMemcpy(d_B, B, size, cudaMemcpyHostToDevice);
36.     VecAdd<<<1, N>>>(d_A, d_B, d_C);
37.     long ed=getCurrentTime();
38.     printf(“gpu running time:%ld\n”,ed–st);
39.     cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);
40.     st=getCurrentTime();
41.     for(int i=0;i<N;i++)
42.         cpu_VecAdd(i,A,B,e);
43.     ed=getCurrentTime();
44.     printf(“cpu running time:%ld\n”,ed–st);
45.     cudaFree(d_A);
46.     cudaFree(d_B);
47.     cudaFree(d_C);
48.  
49.     for(int i=0;i<N;i++)
50.     {
51.         //printf(“%f “,C[i]);
52.     }
53.     printf(“\n”);
54. }

于是

1. make
2. ./test
3. gpu running time:102
4. cpu running time:2263
5. 0.000000 2.000000 8.000000

哇！！GPU好快啊

于是把第十行的j范围增加十倍

1. CUDA$ ./test
2. Hello,World
3. gpu running time:59
4. cpu running time:2121
5. 0.000000 2.000000 8.000000 .

.

是不是毫不科学？

其实是因为，对于内核启动这件事情，实际上是异步的，即我们启动内核以后控制权就交给了host，所以，我们似乎需要一个函数来等待下内核的执行。

其实倒不需要特别修改，因为cudaMemcpy()就是一个同步函数

修改运行代码如下：

1.     cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);
2.     cudaMemcpy(d_B, B, size, cudaMemcpyHostToDevice);
3.     VecAdd<<<1, N>>>(d_A, d_B, d_C);
4.     cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);
5.     long ed=getCurrentTime();
6.     printf(“gpu running time:%ld\n”,ed–st);

执行，将那个毫无意义的j循环设到同等大小

1. Hello,World
2. gpu running time:515
3. cpu running time:2674

又或者提高j循环次数分别提高十倍

1. CUDA$ ./test
2. Hello,World
3. gpu running time:4064
4. cpu running time:24683

这就比较科学啦~

然而需要注意的是内核启动和数据拷贝并非不需要时间，比如如果修改add为如下代码

1. __global__ void VecAdd(float* A, float* B, float* C)
2. {
3.     int i = threadIdx.x;
4.     C[i] = (A[i] * B[i]);
5. }
6. void cpu_VecAdd(int i,float* A, float* B, float* C)
7. {
8.     C[i] = (A[i] * B[i]);
9. }

反倒是GPU慢了

于是，我们已经对CUDA有了初步的尝试

这个Helloworld完整代码如下

1. #include“stdio.h”
2. #include<cuda_runtime.h>
3. #include <sys/time.h>
4.  
5. #define N 1000
6. // Kernel definition
7. __global__ void VecAdd(float* A, float* B, float* C)
8. {
9.     int i = threadIdx.x;
10.     C[i] = (A[i] * B[i]);
11. }
12.  
13. long getCurrentTime()
14. {
15.    struct timeval tv;
16.    gettimeofday(>tv,NULL);
17.    return tv.tv_sec * 1000000 + tv.tv_usec;
18. }
19. void cpu_VecAdd(int i,float* A, float* B, float* C)
20. {
21.     C[i] = (A[i] * B[i]);
22. }
23. int main()
24. {
25.     // Kernel invocation with N threads
26.     printf(“Hello,World\n”);
27.     float *A=new float[N],*B=new float[N],*C=new float[N];
28.     for(int i=0;i<N;i++)
29.     {
30.         A[i]=i;
31.         B[i]=2*i;
32.     }
33.  
34.     size_t size = N * sizeof(float);
35.  
36.     float* d_A;
37.     cudaMalloc(>d_A, size);
38.     float* d_B;
39.     cudaMalloc(>d_B, size);
40.     float* d_C;
41.     cudaMalloc(>d_C, size);
42.     float *e=new float[N];
43.     long st=getCurrentTime();
44.     cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);
45.     cudaMemcpy(d_B, B, size, cudaMemcpyHostToDevice);
46.     VecAdd<<<1, N>>>(d_A, d_B, d_C);
47.     cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);
48.     long ed=getCurrentTime();
49.     printf(“gpu running time:%ld\n”,ed–st);
50.     st=getCurrentTime();
51.     for(int i=0;i<N;i++)
52.         cpu_VecAdd(i,A,B,e);
53.     ed=getCurrentTime();
54.     printf(“cpu running time:%ld\n”,ed–st);
55.     cudaFree(d_A);
56.     cudaFree(d_B);
57.     cudaFree(d_C);
58.  
59.     for(int i=0;i<N;i++)
60.     {
61.         //printf(“%f “,C[i]);
62.     }
63.     printf(“\n”);
64. }