CUDA矩阵元素求和

CUDA矩阵元素求和

---

在使用cuda 做图像高性能处理的时候，会很频繁的使用到矩阵的运算，于是打算做成一个相对完整的矩阵运算库，一旦使用调用即可，便再不用临时去写矩阵运算的规则和方法了。矩阵的加 减 乘 转置自然简单，在cuda中，矩阵的每个元素都对应着一个线程，于是矩阵加法转换为线程中的加法，减法变成了线程，乘法，转置也可以转化成类似的思路。

当我写到求一个矩阵所有元素的和的时候，顺着思维惯性写成了这样的程序

__global__ void matMixKernel(int *Result, int **A){



    //输入参数中，A是要求和的矩阵，Result 是用来保存结果的int指针



    int i = blockIdx.x * blockDim.x + threadIdx.x;

    int j = blockIdx.y * blockDim.y + threadIdx.y;





    Result += A[i][j];





}

抽象一下上面这个程序，仔细一想，自己都笑喷了。

把矩阵想象成一个仓库，大大小小的区域 存放着各种各样的货物，这些货物都用来建房子，

然后每个线程都是一个工人，每个工人管理了一块独立的区域，矩阵的元素求和，就变成了一个很生动的问题
How to build the house？

所以上面的这段代码做得就是：

让所有的工人抱着自己的材料来到要建房子的空地，同时把材料丢下，然后。。。说好的房子就变成了垃圾场，哈哈

不禁想起小学时候老师经常挂在嘴边的问题，十个工人建房子需要一年，四十个工人呢？大家一同喊道，一季度，那么3650个工人呢？又一同喊道，一天。

果真用一天的时间去让3650个人建好一栋房子的话，在五楼干活的那部分人一定会死命的叫唤：老板！四楼水泥还没干呢！

———————————————————分割线———————————————————

在生活中，上面的故事理所应该被当成了笑话，但是在cuda或者是其他种类的并行运算中，类似的笑话经常发生。

把一个有顺序的过程使用并行程序实现，通常被认为是不可能实现的事情，但是稍微转换下思维，并行计算也是可以有顺序存在的。

还是拿刚刚的那个例子，在建房子的过程中，搅和水泥和打地基是可以同时干的两件事情，在这样的子过程中，就是并行计算发挥作用的时候了，我们可以叫100个人同时搅和水泥，叫20人开机器同时打地基，等到地基打好了，水泥也好了，于是乎直接把水泥倒入地基，这样肯定比120个人同时打地基，再搅拌水泥，再倒入水泥速度快。

解决这个问题的方法就是分小队，在cuda中解决类似的问题也是一样，并不是所有线程都需要参与计算，我们把线程分成三种，一种手里拿着加数，一个手里拿着被加数，然后一种手里拿着结果，整个过程也分成两部完成，首先拿着加数的和拿着被加数的碰面做加法运算，算好后手里拿结果的排好队，第二步，把所有手里的结果扔给排头的那个人，然后让那个人拿着总结果来汇报就行了。

是否整个过程行云流水？下面是代码分析，有几个地方我还是改了很久才调试通的

__global__ void block_sum(const int *input,

    int *per_block_results,

    const size_t n)

{

    extern __shared__ int sdata[];



    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;



    // load input into __shared__ memory //一个线程负责把一个元素从全局内存载入到共享内存

    int x = 0;

    if (i < n)

    {

        x = input[i];

    }

    sdata[threadIdx.x] = x;

    __syncthreads();//等待所有线程把自己负责的元素载入到共享内存



    // contiguous range pattern//块内进行合并操作，每次合并变为一半.注意threadIdx.x是块内的偏移，上面算出的i是全局的偏移。

    for (int offset = blockDim.x / 2;

        offset > 0;

        offset >>= 1)

    {

        if (threadIdx.x < offset)//控制只有某些线程才进行操作。

        {

            // add a partial sum upstream to our own

            sdata[threadIdx.x] += sdata[threadIdx.x + offset];

        }



        // wait until all threads in the block have

        // updated their partial sums

        __syncthreads();

    }



    // thread 0 writes the final result//每个块的线程0负责存放块内求和的结果

    if (threadIdx.x == 0)

    {

        per_block_results[blockIdx.x] = sdata[0];

    }

}

//调用方法



cudaError_t matMix(int matA[ROW][COL], unsigned int *Result){



    cudaError_t cudaStatus;



    /***/

    unsigned int num_elements = ROW*COL;

    int h_input[ROW*COL];

    for (int i = 0;i < ROW*COL;i++){



        h_input[i] = matA[i/COL][i%COL];



    }



    //分配内存

    int *d_input = 0;

    cudaMalloc((void**)&d_input, sizeof(int)* num_elements);

    cudaMemcpy(d_input, &h_input[0], sizeof(int)* num_elements, cudaMemcpyHostToDevice);



    const size_t block_size = 512;//线程块的大小。目前有些gpu的线程块最大为512，有些为1024.

    const size_t num_blocks = (num_elements / block_size) + ((num_elements%block_size) ? 1 : 0);





    int *d_partial_sums_and_total = 0;//一个线程块一个和，另外加一个元素，存放所有线程块的和。

    cudaMalloc((void**)&d_partial_sums_and_total, sizeof(int)* (num_blocks + 1));



//把每个线程块的和求出来

    block_sum <<<num_blocks, block_size, block_size *sizeof(int)>>>(d_input, d_partial_sums_and_total, num_elements);



    cudaStatus = cudaThreadSynchronize();





    //再次用一个线程块把上一步的结果求和。

    //注意这里有个限制，上一步线程块的数量，必须不大于一个线程块线程的最大数量，因为这一步得把上一步的结果放在一个线程块操作。

    //即num_blocks不能大于线程块的最大线程数量。

    block_sum << <1, num_blocks, num_blocks * sizeof(int) >> >(d_partial_sums_and_total, d_partial_sums_and_total + num_blocks, num_blocks);

    cudaStatus = cudaThreadSynchronize();





    // copy the result back to the host

    int device_result = 0;

    cudaMemcpy(&device_result, d_partial_sums_and_total + num_blocks, sizeof(int), cudaMemcpyDeviceToHost);



    *Result = device_result;



    if (cudaStatus != cudaSuccess){



        cout <<cudaStatus << endl;



    }



    // deallocate device memory

    cudaFree(d_input);

    cudaFree(d_partial_sums_and_total);









    return cudaStatus;

}

于是调试，搞定！