[CUDA]CUDA编程实战五——dot向量点积运算

在前面是有考虑把点积运算放到前面，但是考虑到有一些新的东西，还是准备把他放在后面了。

问题说明

对于向量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，和前一个优化方案的时间大致相同。