在前面是有考虑把点积运算放到前面，但是考虑到有一些新的东西，还是准备把他放在后面了。
问题说明
对于向量a[1-n]以及b[1-n]，点积运算即为将对应位置上的元素相加，即c = a[1]b[1]+...a[i]b[i]+...+a[n]b[n];
为了方便起见，我们将设置a[i]=i,b[i]=2*i;那么其最终结果为平方和的两倍。
简单CUDA版
对于长度为N的vector，我们可以开一个block，每个block中有1024个线程，这样每个线程都将执行一次乘加运算，最后再将结果累加。
这里我们使用了共享内存，共享内存是同一个block中概念，同一个block中的线程可以通过共享内存的方式进行线程间通信。
每个线程计算一个point，比矩阵乘法那种每个线程计算一行的方式要快很多，最后使用线程0进行收缩。
#include "cuda_runtime.h"
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
#include <cstdio>
using namespace std;
const int N = 1024;
float sum_squares(float x)
{
return x * (x+1) * (2 * x + 1) / 6;
}
__global__
void dot_simple(float* a, float* b, float* c)
{
__shared__ float temp[N];
int i = threadIdx.x + blockIdx.x * blockDim.x;
temp[i] = a[i] * b[i];
__syncthreads();
if (0 == i)
{
int sum = 0;
for (int i=0; i < N; i++)
{
sum += temp[i];
}
*c = sum;
printf("sum Calculated on Device: %.6g
", *c);
}
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
float* a;
float* b;
float* c;
float* A;
float* B;
float* C;
int sz = N * sizeof(float);
a = (float*)malloc(sz);
b = (float*)malloc(sz);
c = (float*)malloc(sz);
cudaMalloc(&A, sz);
cudaMalloc(&B, sz);
cudaMalloc(&C, sz);
for (int i=0; i<N; i++)
{
a[i] = i;
b[i] = i * 2;
}
cudaMemcpy(A, a, sz, cudaMemcpyHostToDevice);
cudaMemcpy(B, b, sz, cudaMemcpyHostToDevice);
dim3 blocks(1);
dim3 threads(1024);
dot_simple<<<blocks, threads>>>(A, B, C);
cudaMemcpy(c, C, sz, cudaMemcpyDeviceToHost);
float sum_c = c[0];
// for (int i=0; i<N/256; i++)
// {
// sum_c += c[i];
// }
printf("GPU results: %.6g == %.6g
", sum_c, 2*sum_squares((float)N-1));
free(a);
free(b);
free(c);
cudaFree(A);
cudaFree(B);
cudaFree(C);
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
}
运行结果
我们将cuda运行结果和cpu的运行结果进行对比，结果正确，最快花费了108ms。
cuda优化版
上述的实现存在一定问题，即只能在一个block中进行，但每个block最大线程数是1024，那么对应向量的长度也最长为1024，极大约束了算法的扩展。
因此我们做了一些优化，使用256个线程，每个线程的计算结果是i+256*k的累积结果，然后使用cpu将这些结果再累加起来。
这里我们使用了一个新的知识点，原子累加，即通过原子操作来防止加的过程中不同步问题,注意我们使用的是atomicadd，而不是+="，因为后者会带来不同步的问题。
__global__
void dot2(float* a, float* b, float* c)
{
int i = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIdx = threadIdx.x;
float sum = 0.0;
while (i < N){
sum += a[i] * b[i];
i += blockDim.x * gridDim.x;
}
// c[cacheIdx] += sum; //会出现不同步的问题
atomicAdd(&c[cacheIdx], sum);
__syncthreads();
// printf("%d %d %d %d
", i, threadIdx.x, blockIdx.x, blockDim.x);
}
dim3 blocks(N/256);
dim3 threads(256);
dot2<<<blocks, threads>>>(A, B, C);
cudaMemcpy(c, C, sz, cudaMemcpyDeviceToHost);
float sum_c = 0.0;
for (int i=0; i<256; i++)
{
sum_c += c[i];
}
运行结果
我们发现运行结果有了显著的提升，在最好的情况下，运行时间为90ms。
CUDA优化版2
在上一个版本中，将加法运算放到外面，有没有更好的操作，显然有的，快速的方案是将进行规约的过程也进行多线程化，多个线程进行最后的加法，需要logn次，而非n次。
需要注意的是， 第二个__syncthreads()必须写在if外面，因为 __syncthreads()是针对所有线程的，如果写在里面会引发线程阻塞并使得程序崩溃。
__global__
void dot3(float* a, float* b, float* c)
{
const int threadsPerBlock = 256;
__shared__ float cache[threadsPerBlock];
int i = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIdx = threadIdx.x;
float sum = 0;
while (i < N){
sum += a[i] * b[i];
i += blockDim.x * gridDim.x;
}
cache[cacheIdx] = sum;
__syncthreads();
int num = blockDim.x / 2;
while (num!=0){
if (cacheIdx < num) cache[cacheIdx] = cache[cacheIdx] + cache[cacheIdx + num];
__syncthreads();
num /= 2;
}
if (cacheIdx == 0) c[blockIdx.x] = cache[0];
}
dim3 blocks(N/256);
dim3 threads(256);
dot2<<<blocks, threads>>>(A, B, C);
cudaMemcpy(c, C, sz, cudaMemcpyDeviceToHost);
float sum_c = 0.0;
for (int i=0; i<N/256; i++)
{
sum_c += c[i];
}
运行结果
最快的运行结果为92ms，和前一个优化方案的时间大致相同。