在使用cuda 作图像高性能处理的时候，会很频繁的使用到矩阵的运算，因而打算作成一个相对完整的矩阵运算库，一旦使用调用便可，便再不用临时去写矩阵运算的规则和方法了。矩阵的加 减 乘 转置天然简单，在cuda中，矩阵的每一个元素都对应着一个线程，因而矩阵加法转换为线程中的加法，减法变成了线程，乘法，转置也能够转化成相似的思路。web
当我写到求一个矩阵全部元素的和的时候，顺着思惟惯性写成了这样的程序svg
__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];
}
抽象一下上面这个程序，仔细一想，本身都笑喷了。
性能
把矩阵想象成一个仓库，大大小小的区域 存放着各类各样的货物，这些货物都用来建房子，ui
而后每一个线程都是一个工人，每一个工人管理了一块独立的区域，矩阵的元素求和，就变成了一个很生动的问题
How to build the house？
线程
因此上面的这段代码作得就是：
指针
让全部的工人抱着本身的材料来到要建房子的空地，同时把材料丢下，而后。。。说好的房子就变成了垃圾场，哈哈
调试
不由想起小学时候老师常常挂在嘴边的问题，十个工人建房子须要一年，四十个工人呢？你们一同喊道，一季度，那么3650个工人呢？又一同喊道，一天。
code
果然用一天的时间去让3650我的建好一栋房子的话，在五楼干活的那部分人必定会死命的叫唤：老板！四楼水泥还没干呢！
xml
———————————————————分割线———————————————————
内存
在生活中，上面的故事理所应该被当成了笑话，可是在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;
}
因而调试，搞定！