在前面是有考慮把點積運算放到前面，但是考慮到有一些新的東西，還是準備把他放在後面了。

問題說明

對於向量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\n", *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\n", 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\n", 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，和前一個優化方案的時間大致相同。