问题：
给定两个尺寸为N的向量a和b，
其中对于从0到N的每一个i，a[i] = i，b[i] = i，
计算二者点积。
参考结果：
首先给出python中numpy的计算结果，为了方便比较，我输出了用int类型和float类型计算的两种结果。
当N比较小时，以N = 100为例：
$ python dot_product.py
328350.0
328350
当N稍微大一点，以N = 1,000为例：
$ python dot_product.py
332833500.0
332833500
当N = 1,000,000时：
$ python dot_product.py
3.3333283333349997e+17
333332833333500000
可以看到即使是numpy这样十分完善的科学计算库，也会因为精度不够而导致结果有偏差的问题。
python程序源代码
import numpy as np
N = 1000000
a = []
b = []
c = []
d = []
for i in range(N):
a.append(float(i))
b.append(float(i))
a = np.array(a)
b = np.array(b)
e = np.dot(a,b)
print(e)
for i in range(N):
c.append(i)
d.append(i)
c = np.array(c)
d = np.array(d)
f = np.dot(c,d)
print(f)
C语言下CPU和GPU计算结果比较：
首先来看以float类型来计算
N = 100：
$ nvcc -o lab02_2_my.o lab02_2_my.cu
$ ./lab02_2_my.o
GPU: 328350.000000
CPU: 328350.000000
结果与之前numpy计算结果一致。
N = 1,000：
$ nvcc -o lab02_2_my.o lab02_2_my.cu
$ ./lab02_2_my.o
GPU: 332833408.000000
CPU: 332833152.000000
这时的结果相比numpy已经有些偏差了，而且是GPU和CPU的计算结果也出现了不同。
N = 1,000,000：
$ nvcc -o lab02_2_my.o lab02_2_my.cu
$ ./lab02_2_my.o
GPU: 333379003647787008.000000
CPU: 333380996512612352.000000
可以看到此时的偏差越来越大。
再来看以int类型来计算
N = 100：
$ nvcc -o lab02_2_my_int.o lab02_2_my_int.cu
$ ./lab02_2_my_int.o
GPU: 328350
CPU: 328350
结果与之前numpy计算结果一致。
N = 1,000：
$ nvcc -o lab02_2_my_int.o lab02_2_my_int.cu
$ ./lab02_2_my_int.o
GPU: 332833500
CPU: 332833500
这时的结果同样与之前numpy相同，且在CPU和GPU上计算的结果也相同。
N = 1,000,000：
$ nvcc -o lab02_2_my_int.o lab02_2_my_int.cu
$ ./lab02_2_my_int.o
GPU: 584144992
CPU: 584144992
这个结果就差得远了，连数量级都不对。
补充一个 N = 10,000的结果
$ nvcc -o lab02_2_my_int.o lab02_2_my_int.cu
$ ./lab02_2_my_int.o
GPU: -1724114088
CPU: -1724114088
从这个结果上看，应该是数值溢出了。
但是！请注意：虽然结果不对，但是这次CPU和GPU的运算结果是一致的！
重新检查了一下代码，发现我为了debug把GPU的block数和thread数都设置成了1，也就是说在GPU里只开辟了一个thread来进行计算，这是唯一与之前float程序的不同之处。为了验证我把float程序的也修改成了block = 1，thread = 1， 结果如下：
$ nvcc -o lab02_2_my.o lab02_2_my.cu
$ ./lab02_2_my.o
GPU: 333380996512612352.000000
CPU: 333380996512612352.000000
这次CPU和GPU的结果也一致了！然后我又试了以下几种组合：
Block number
Thread number
GPU == CPU?
1
1
true
1
64
true
1
128
false
2
1
true
2
32
true
2
64
false
4
16
true
4
32
false
16
1
true
16
2
false
从以上结果可以看出，block和thread的数量影响着GPU计算的精度，而二者相比block的影响更为显著，我认为结果不一致可能是显存相关的原因所致，具体原因我需要再问一问professor。
又试了一下设置block和thread为dim3类型的参数，又出现无法理解的结果了：
const dim3 BLOCKS_NUM = dim3(1,1,1);
const dim3 THREADS_NUM = dim3(2,2,1);
结果：
$ nvcc -o lab02_2_my.o lab02_2_my.cu
$ ./lab02_2_my.o
GPU: 666259653650284544.000000
CPU: 333380996512612352.000000
出现了将近2倍的偏差，思考了一下好像是核函数中Grid-Striding Loops那部分只考虑到了x轴，要想引入y轴和z轴，需要将代码修改为
for (i = blockIdx.x * blockDim.x + blockIdx.y * blockDim.y + blockIdx.z * blockDim.z+ threadIdx.x + threadIdx.y + threadIdx.z;i < N;i += blockDim.x * gridDim.x * blockDim.y * gridDim.y*blockDim.z * gridDim.z)
{
tmp_product = a[i] * b[i];
atomicAdd(c, tmp_product);
}
结果也就正常了：
$ nvcc -o lab02_2_my.o lab02_2_my.cu
$ ./lab02_2_my.o
GPU: 333357494451568640.000000
CPU: 333380996512612352.000000
结论：
GPU对大数的运算，受所开辟的block和thread数量影响，当block与thread少时，无法体现出GPU的运算优势，而数量多时，又会影响计算结果的精度。
C程序源代码
//hubutang 15.10.2017 点积计算：float 类型
#include
#define CUDA_CHECK(cmd) {cudaError_t error = cmd; if(error!=cudaSuccess){printf(":%i ",__FILE__,__LINE__); printf("[CUDA] Error: %s\n", cudaGetErrorString(error));}}
#define CUDA_CHECK_KERNEL {cudaError_t error = cudaGetLastError(); if(error!=cudaSuccess){printf(":%i ",__FILE__,__LINE__); printf("[CUDA] Error: %s\n", cudaGetErrorString(error));}}
//以上两行是为了检测GPU操作中可能出现的错误
#define N 10000
//数组的尺寸
//======================== part1. GPU计算点积的核函数 ======================//
__global__ void kernel_dot_product( float *a, float *b, float *c )
{
int i;
float tmp_product;
for (i = blockIdx.x * blockDim.x + threadIdx.x;i < N;i += blockDim.x * gridDim.x)
//这个for循环是为了提高效率：Grid-Striding Loops
{
tmp_product = a[i] * b[i];
atomicAdd(c, tmp_product);
//cuda中的原子操作，涉及到需要对共享内存中
//的内容操作时使用，本例中既是求和的步骤
}
}
//========================================================================//
//======================== part2. CPU计算点积的核函数 ======================//
float cpu_dot_product( float *a, float *b)
{
int i=0;
float c=0.0;
while(1)
{
c = c + a[i] * b[i];
i++;
if( i>=N )
{
break;
}
}
return c;
}
//========================================================================//
//===================== part3. 主函数 比较两个函数的输出结果 ==================//
int main (int argc, char **argv){
float *a_host = new float[N];
float *b_host = new float[N];
float *c_host = new float;
float zero = 0.0;
c_host = &zero;
float *a_device, *b_device, *c_device;
int i;
for (i=0;i
{
a_host[i]=float(i);
b_host[i]=float(i);
}
CUDA_CHECK(cudaMalloc(&a_device, N*sizeof(float)));
CUDA_CHECK(cudaMalloc(&b_device, N*sizeof(float)));
CUDA_CHECK(cudaMalloc(&c_device, 1*sizeof(float)));
//在GPU中分配内存
CUDA_CHECK(cudaMemcpy(a_device, a_host, N*sizeof(float), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(b_device, b_host, N*sizeof(float), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(c_device, c_host, 1*sizeof(float), cudaMemcpyHostToDevice));
//将数据传输到GPU上
const int BLOCKS_NUM = 100;
const int THREADS_NUM = 512;
//GPU参数设置：使用100个block，每个block开启512个thread
kernel_dot_product<<>>(a_device, b_device, c_device);
//启动GPU核函数，计算点积
CUDA_CHECK(cudaMemcpy(c_host, c_device, 1*sizeof(float), cudaMemcpyDeviceToHost));
//计算结果传输到主机，以便打印输出
printf("GPU: %f \n",*c_host);
float cpu_result = cpu_dot_product(a_host, b_host);
//CPU计算点积
printf("CPU: %f \n",cpu_result);
CUDA_CHECK(cudaFree(a_device));
CUDA_CHECK(cudaFree(b_device));
CUDA_CHECK(cudaFree(c_device));
CUDA_CHECK( cudaDeviceReset() );
delete[] a_host;
delete[] b_host;
//释放内存
return 0;
}
//========================================================================//
//hubutang 15.10.2017 点积计算：int类型
#include
#define CUDA_CHECK(cmd) {cudaError_t error = cmd; if(error!=cudaSuccess){printf(":%i ",__FILE__,__LINE__); printf("[CUDA] Error: %s\n", cudaGetErrorString(error));}}
#define CUDA_CHECK_KERNEL {cudaError_t error = cudaGetLastError(); if(error!=cudaSuccess){printf(":%i ",__FILE__,__LINE__); printf("[CUDA] Error: %s\n", cudaGetErrorString(error));}}
//以上两行是为了检测GPU操作中可能出现的错误
#define N 10000
//数组的尺寸
//======================== part1. GPU计算点积的核函数 ======================//
__global__ void kernel_dot_product( int *a, int *b, int *c )
{
int i;
int tmp_product;
for (i = blockIdx.x * blockDim.x + threadIdx.x;i < N;i += blockDim.x * gridDim.x)
//这个for循环是为了提高效率：Grid-Striding Loops
{
tmp_product = a[i] * b[i];
atomicAdd(c, tmp_product);
//cuda中的原子操作，涉及到需要对共享内存中
//的内容操作时使用，本例中既是求和的步骤
}
}
//========================================================================//
//======================== part2. CPU计算点积的核函数 ======================//
int cpu_dot_product( int *a, int *b)
{
int i=0;
int c=0;
while(1)
{
c = c + a[i] * b[i];
i++;
if( i>=N )
{
break;
}
}
return c;
}
//========================================================================//
//===================== part3. 主函数 比较两个函数的输出结果 ==================//
int main (int argc, char **argv){
int *a_host = new int[N];
int *b_host = new int[N];
int *c_host = new int;
int zero = 0;
c_host = &zero;
int *a_device, *b_device, *c_device;
int i;
for (i=0;i
{
a_host[i]=i;
b_host[i]=i;
}
CUDA_CHECK(cudaMalloc(&a_device, N*sizeof(int)));
CUDA_CHECK(cudaMalloc(&b_device, N*sizeof(int)));
CUDA_CHECK(cudaMalloc(&c_device, 1*sizeof(int)));
//在GPU中分配内存
CUDA_CHECK(cudaMemcpy(a_device, a_host, N*sizeof(int), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(b_device, b_host, N*sizeof(int), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(c_device, c_host, 1*sizeof(int), cudaMemcpyHostToDevice));
//将数据传输到GPU上
const int BLOCKS_NUM = 1;
const int THREADS_NUM = 1;
//GPU参数设置：使用1个block，每个block开启1个thread
kernel_dot_product<<>>(a_device, b_device, c_device);
//启动GPU核函数，计算点积
CUDA_CHECK(cudaMemcpy(c_host, c_device, 1*sizeof(int), cudaMemcpyDeviceToHost));
printf("GPU: %d \n",*c_host);
int cpu_result = cpu_dot_product(a_host, b_host);
//CPU计算点积
printf("CPU: %d \n",cpu_result);
CUDA_CHECK(cudaFree(a_device));
CUDA_CHECK(cudaFree(b_device));
CUDA_CHECK(cudaFree(c_device));
CUDA_CHECK( cudaDeviceReset() );
delete[] a_host;
delete[] b_host;
//释放内存
return 0;
}
//========================================================================//
