问题：

给定两个尺寸为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和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 <stdio.h>

#define CUDA_CHECK(cmd) {cudaError_t error = cmd; if(error!=cudaSuccess){printf("<%s>:%i ",__FILE__,__LINE__); printf("[CUDA] Error: %s\n", cudaGetErrorString(error));}}
#define CUDA_CHECK_KERNEL {cudaError_t error = cudaGetLastError(); if(error!=cudaSuccess){printf("<%s>:%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<N;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<<<BLOCKS_NUM, THREADS_NUM>>>(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 <stdio.h>

#define CUDA_CHECK(cmd) {cudaError_t error = cmd; if(error!=cudaSuccess){printf("<%s>:%i ",__FILE__,__LINE__); printf("[CUDA] Error: %s\n", cudaGetErrorString(error));}}
#define CUDA_CHECK_KERNEL {cudaError_t error = cudaGetLastError(); if(error!=cudaSuccess){printf("<%s>:%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<N;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<<<BLOCKS_NUM, THREADS_NUM>>>(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;
}
//========================================================================//