CUDA学习笔记三_牛客博客

我的代码：https://github.com/BBuf/CUDA_LEARN_SAMPLES/tree/58b95d34b60c9b50c72d66a2d2b3b671337191bd
CUDA加速的非常核心的一个应用是对矩阵运算进行加速，所以这篇博客，我会记录利用CUDA写矩阵乘法加速的过程。
矩阵乘法 所谓矩阵乘法就是了AB=C，其中A为n*m的矩阵,B为m*q的矩阵,C为n*q的矩阵，矩阵乘法的代码如下：

const LL mod = 1e9+7;
struct Matrix{
    LL a[2][2];
    void set1(){
        memset(a, 0, sizeof(a));
    }
    void set2(){
        set1();
        for(int i=0; i<2; i++) a[i][i]=1;
    }
};
Matrix operator*(const Matrix &a, const Matrix &b){
    Matrix res;
    res.set1();
    for(int i=0; i<2; i++){
        for(int j=0; j<2; j++){
            for(int k=0; k<2; k++){
                res.a[i][j] = (res.a[i][j] + a.a[i][k]*b.a[k][j]%mod)%mod;
            }
        }
    }
    return res;
}

这里为了方便，接下来涉及的矩阵均为n*n的方形矩阵。
并行矩阵乘法核函数

//n * n Matrix Multiply A * B = C
__global__ static void MatrixMulCUDA(const float *a, const float *b, float *c, int n)
{
    const int tid = threadIdx.x; //目前的thread是第几个thread
    const int bid = blockIdx.x; //目前的thread是属于哪个block
    const int idx = bid * THREAD_NUM + tid; //从bid和tid推出目前的thread应该计算的A矩阵的行数和B矩阵列数
    const int row = idx / n;
    const int col = idx % n;
    if(row < n && col < n){
        float sum = 0;
        for(int i = 0; i < n; i++){
            sum += a[row * n + i] * b[i * n + col];
        }
        c[row * n + col] = sum;
    }
}

生成随机矩阵

//生成随机矩阵
void RandomMatrix(float *A, int n){
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            A[i * n + j] = (float)rand() / MY_RAND_MAX + (float)rand() / (MY_RAND_MAX * MY_RAND_MAX);
        }
    }
}

矩阵运算

void MatrixMul(){
    //定义矩阵
    float *a, *b, *c, *d;
    int n = MATRIX_SIZE;
    a = (float*)malloc(sizeof(float)*n*n);
    b = (float*)malloc(sizeof(float)*n*n);
    c = (float*)malloc(sizeof(float)*n*n);
    d = (float*)malloc(sizeof(float)*n*n);
    srand(time(NULL));
    RandomMatrix(a, n);
    RandomMatrix(b, n);
    //把数据复制到显卡内存中
    float *cuda_a, *cuda_b, *cuda_c;
    //cudaMalloc 取得一块显存内存
    cudaMalloc((void**)&cuda_a, sizeof(float) * n * n);
    cudaMalloc((void**)&cuda_b, sizeof(float) * n * n);
    cudaMalloc((void**)&cuda_c, sizeof(float) * n * n);
    //cudaMemcpy 将产生的矩阵复制到显卡内存中
    //cudaMemcpyHostToDevice - 从内存复制到显卡内存
    //cudaMemcpyDeviceToHost - 从显卡内存复制到内存
    cudaMemcpy(cuda_a, a, sizeof(float)*n*n, cudaMemcpyHostToDevice);
    cudaMemcpy(cuda_b, b, sizeof(float)*n*n, cudaMemcpyHostToDevice);
    float time_elapsed=0;
    cudaEvent_t start,stop;
    cudaEventCreate(&start);    //创建Event
    cudaEventCreate(&stop);
    cudaEventRecord( start,0);    //记录当前时间
    MatrixMulCUDA<<<blocks_num, THREAD_NUM, 0>>>(cuda_a, cuda_b, cuda_c, n);
    cudaEventRecord( stop,0);    //记录当前时间
    cudaEventSynchronize(start);    //Waits for an event to complete.
    cudaEventSynchronize(stop);    //Waits for an event to complete.Record之前的任务
    //把结果从显示芯片复制回主内存
    //cudaMemcpy 将结果从显存中复制回内存
    cudaMemcpy(c, cuda_c, sizeof(float) * n * n, cudaMemcpyDeviceToHost);
    cudaEventElapsedTime(&time_elapsed,start,stop);    //计算时间差
    cudaEventDestroy(start);    //destory the event
    cudaEventDestroy(stop);
    //Free
    cudaFree(cuda_a);
    cudaFree(cuda_b);
    cudaFree(cuda_c);
    printf("Matrix multiply GPU time: %.10f\n", time_elapsed);
    clock_t start_time = clock();
    //CPU计算矩阵乘法
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            double temp = 0;
            for(int k = 0; k < n; k++){
                temp += a[i * n + k]* b[k * n + j];
            }
            d[i * n + j] = temp;
        }
    }
    clock_t end_time = clock();
    //验证正确性和准确性
    float max_error = 0.0, average_error = 0;
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            if(d[i * n + j] != 0){
                float err = fabs((c[i * n + j] - d[i * n + j]) / d[i * n + j]);
                if(max_error < err) max_error = err;
                average_error += err;
            }
        }
    }
    double cpu_time = (double)(end_time - start_time) / CLOCKS_PER_SEC * 1000.0;
    printf("Matrix multiply CPU time: %.10f\n", cpu_time);
    printf("Max error: %.10f Average error: %.10f\n", max_error, average_error / (n * n));
    printf("%.10f\n", cpu_time/time_elapsed);
}

矩阵运算出现精度问题 由于在 CPU 上进行计算时，我们使用 double（即 64 bits 浮点数）来累进计算过程，而在 GPU 上则只能用 float（32 bits 浮点数）。在累加大量数字的时候，由于累加结果很快会变大，因此后面的数字很容易被舍去过多的位数。例如下面这段代码：

    float a = 100998;
    float b = 2.338;
    a = a + b;
    printf("the sum is %.10f\n", a);

输出是：101000.3359375000，容易看出有一个0.008被截断了，在1000次操作后一个整数8就没了，所以造成了很大的精度问题。
Kahan’s Summation Formula 为了一定程度上解决这个精度问题，我们需要记住这个损失的误差值，假设误差为temp= (a+b)-a-b, 在上面那个问题中 temp = -0.008,在下次计算的时候加和到下一个加数就可以一定程度的减小误差。
有一篇这个算法的好文章：https://www.cnblogs.com/kid551/p/6444520.html 这个算法的伪代码描述如下：

def KahanSum(input):
    var sum = 0.0
    var c = 0.0
    for i = 1 to input.length do
        var y = input[i] - c    // Initially, c is zero; then it compensates previous accuracy.
        var t = sum + y         // low-order digits of y are lost
        c = (t - sum) - y       // recover the low-order digits of y, with negative symbol
        sum = t
    next i

    return sum

在上述伪代码中，变量c表示的即是小数的补全部分compensation，更严格地说，应该是负的补全部分。随着这个补全部分的不断积累，当这些截断误差积累到一定量级，它们在求和的时候也就不会被截断了，从而能够相对好地控制整个求和过程的精度。然后把这个算法应用到我们的矩阵乘法核函数中。

//n * n Matrix Multiply A * B = C
__global__ static void MatrixMulCUDA(const float *a, const float *b, float *c, int n)
{
    const int tid = threadIdx.x; //目前的thread是第几个thread
    const int bid = blockIdx.x; //目前的thread是属于哪个block
    const int idx = bid * THREAD_NUM + tid; //从bid和tid推出目前的thread应该计算的A矩阵的行数和B矩阵列数
    const int row = idx / n;
    const int col = idx % n;
    if(row < n && col < n){
// float sum = 0;
// for(int i = 0; i < n; i++){
// sum += a[row * n + i] * b[i * n + col];
// }
        float t = 0;
        float y = 0;
        for(int i = 0; i < n; i++){
            float r;
            y -= a[row * n + i] * b[i * n + col];
            r = t - y;
            y = (r - t) + y;
            t = r;
        }
        c[row * n + col] = t;
    }
}

然后不知道为什么，在我的GPU上的精度是没有提升的。。。