一、Conv
实现一个卷积操作，使用c语言编写的话实际就是四个循环，一个循环遍历输入矩阵的行，一个循环遍历输入矩阵的列，定位到需要卷积的子矩阵，然后一个循环遍历kernel的行，一个循环遍历kernel的列，这两个循环用来计算当前卷积的结果。
1.基础版本
void conv2d(float *input,float *kernel,float *output,const int in_size,const int ker_size,const int stride){
int out_size=(in_size-ker_size)/stride+1;
for(int i=0;i<out_size;++i){
for(int j=0;j<out_size;++j){
float sum=0;
for(int ker_i=0;ker_i<ker_size;++ker_i){
for(int ker_j=0;ker_j<ker_size;++ker_j){
sum+=kernel[ker_i*ker_size+ker_j]*input[(i*stride+ker_i)*in_size+(j*stride+ker_j)];
}
}
output[i*out_size+j]=sum;
}
}
}
2.neon版本
对于直接使用neon函数，以及直接使用滑动窗口的方式进行卷积的话，因为每一个neon寄存器存储4个float数据，并且深度学习中使用的大多数都是奇数行列的卷积核（比如3*3），所以直接使用滑动窗口的方式进行卷积会浪费掉neon的空间。为了能够让neon寄存器充分的利用起来，可以使用im2col的方式将需要卷积的子矩阵统一变层一维的一个列，这样就可以使用之前实现的矩阵乘来进行计算了。
关于im2col的原理，可以参考其他博客。
总体来说，im2col的意思就是说将输入矩阵按照需要卷积的对应子矩阵转换为一维的向量，很多地方的转换是将输入矩阵的子矩阵转换为一个列，本实现转换为一个行。本实验使用的输入矩阵是640*640，卷积核大小为3*3，步长为1。所以按照im2col转换的话，输入会转换为一个638*638行、3*3列的一个矩阵。im2col代码如下：
void im2col(const float *input,float *output,const int in_size,const int ker_size,const int h,const int w,const int stride){
int out_size=(in_size-ker_size)/stride+1;
for(int i=0;i<out_size;++i){
for(int j=0;j<out_size;++j){
for(int ker_i=0;ker_i<ker_size;++ker_i){
for(int ker_j=0;ker_j<ker_size;++ker_j){
int in_index=(i*stride+ker_i)*in_size+j*stride+ker_j;
int out_index=(i*out_size+j)*ker_size*ker_size+ker_size*ker_i+ker_j;
output[out_index]=input[in_index];
}
}
}
}
}
在将输入矩阵转换后，就可以使用neon intrinsic函数来并行处理数据了。由于我们将每一个需要卷积的矩阵都转换为了一个行，那进行一个卷积操作就是对一行进行一个乘加操作，使用卷积核对整个转换后的矩阵的每一行都进行乘加操作后，就完成了对整张图像的卷积操作了。
下面是整体代码，仅供参考：
void conv2d(const float *input,const float *kernel,const int in_size,const int ker_size,const int stride,float *output){
int out_size=(in_size-ker_size)/stride+1;
int w=ker_size*ker_size;
int h=out_size*out_size;
float *input_col=(float*)malloc(sizeof(float)*h*w);
im2col(input,input_col,in_size,ker_size,h,w,stride);
int w_len=(w/4)*4;
int h_len=(h/4)*4;
for(int i=0;i<h;i+=4){
float32x4_t sum0=vdupq_n_f32(0);
float32x4_t sum1=vdupq_n_f32(0);
float32x4_t sum2=vdupq_n_f32(0);
float32x4_t sum3=vdupq_n_f32(0);
float res0=0;
float res1=0;
float res2=0;
float res3=0;
int j;
for(j=0;j<w_len;j+=4){
float32x4_t va=vld1q_f32(&kernel[j]);
float32x4_t vb0=vld1q_f32(&input_col[i*w+j]);
float32x4_t vb1=vld1q_f32(&input_col[(i+1)*w+j]);
float32x4_t vb2=vld1q_f32(&input_col[(i+2)*w+j]);
float32x4_t vb3=vld1q_f32(&input_col[(i+3)*w+j]);
sum0=vmlaq_f32(sum0,va,vb0);
sum1=vmlaq_f32(sum1,va,vb1);
sum2=vmlaq_f32(sum2,va,vb2);
sum3=vmlaq_f32(sum3,va,vb3);
}
float32x2_t sum_part0=vadd_f32(vget_low_f32(sum0),vget_high_f32(sum0));
sum_part0=vpadd_f32(sum_part0,sum_part0);
res0=vget_lane_f32(sum_part0,0);
float32x2_t sum_part1=vadd_f32(vget_low_f32(sum1),vget_high_f32(sum1));
sum_part1=vpadd_f32(sum_part1,sum_part1);
res1=vget_lane_f32(sum_part1,0);
float32x2_t sum_part2=vadd_f32(vget_low_f32(sum2),vget_high_f32(sum2));
sum_part2=vpadd_f32(sum_part2,sum_part2);
res2=vget_lane_f32(sum_part2,0);
float32x2_t sum_part3=vadd_f32(vget_low_f32(sum3),vget_high_f32(sum3));
sum_part3=vpadd_f32(sum_part3,sum_part3);
res3=vget_lane_f32(sum_part3,0);
if(j<w){
for(;j<w;++j){
res0+=input_col[i*w+j]*kernel[j];
res1+=input_col[(i+1)*w+j]*kernel[j];
res2+=input_col[(i+2)*w+j]*kernel[j];
res3+=input_col[(i+3)*w+j]*kernel[j];
}
}
output[i]=res0;
output[i+1]=res1;
output[i+2]=res2;
output[i+3]=res3;
}
free(input_col);
}
在不打开O3优化的运行时间如下：
在打开03优化后的运行时间为：
3.cuda版本
使用cuda处理卷积时，按照每一个线程处理一个子卷积处理，具体代码如下：
__global__
void conv2d(const float *input,const float *kernel,const int in_size,const int ker_size,const int stride,float *output){
int row=blockIdx.y*blockDim.y+threadIdx.y;
int col=blockIdx.x*blockDim.x+threadIdx.x;
const int out_size=(in_size-ker_size)/stride+1;
float sum=0;
if(row<out_size&&col<out_size){
for(int i=0;i<ker_size;++i){
for(int j=0;j<ker_size;++j){
sum+=input[(row+i)*in_size+(col+j)]*kernel[i*ker_size+j];
}
}
output[row*out_size+col]=sum;
}
__syncthreads();
}
使用nvprof分析如下：
使用cuda处理输入为640*640,、卷积核大小为3*3，步长为1时消耗时间为13.171ms，即使加上数据传输也才不到17ms，与使用neon intrinsic函数处理快了10ms左右。
二、Pool
1.maxpool
（1）基础版本
maxpool的计算原理很简单，对于输入张量，根据给出的kernel大小，选择出输入张量的对应kernel内部的最大值作为当前kernel计算的结果。
以下是使用c++实现的基础版本。一个循环遍历输入数据的行，一个循环遍历当前行的列，然后针对定位到的每一个行和每一个列都遍历该对应kernel内的数据，找出该kernel中的最大值作为当前kernel计算的结果。具体代码如下：
void naive(float *input,const int in_size,const int ker_size,const int stride,float *output){
int out_size=(in_size-ker_size)/stride+1;
for(int i=0;i<out_size;++i){
for(int j=0;j<out_size;++j){
int max=0;
for(int ker_i=0;ker_i<ker_size;++ker_i){
for(int ker_j=0;ker_j<ker_size;++ker_j){
int row=i*stride+ker_i;
int col=j*stride+ker_j;
max=max<input[row*in_size+col]?input[row*in_size+col]:max;
}
}
output[i*out_size+j]=max;
}
}
}
输入数据为640*640的矩阵，池化核大小为3*3，步长为1时，不打开O3运行时间如下：
打开O3优化后的运行时间为：
上述实现中，不打开O3优化时运行时间为94ms，打开O3优化运行时间为29ms左右。
（2）neon版本
对于每一个需要池化的核大小的矩阵，在选取最大值的时候可以使用neon intrinsic函数进行优化。对于3*3的池化操作，具体来说可以使用三个寄存器对每一个池化的矩阵进行存储，然后使用vdupq_n_f32对每一个寄存器的第4个数据赋值为0，那么就可以保证操作的数据位3*3的矩阵。
然后对这三个寄存器进行比较，找出对应位置的最大值，最后获取结果寄存器内部元素的最大值。
具体代码如下：
void maxpool(float *input,const int in_size,const int ker_size,const int stride,float *output){
int out_size=(in_size-ker_size)/stride+1;
for(int i=0;i<out_size;++i){
for(int j=0;j<out_size;++j){
float32x4_t va=vdupq_n_f32(0);
float32x4_t vb=vdupq_n_f32(0);
float32x4_t vc=vdupq_n_f32(0);
float32x4_t max_v=vdupq_n_f32(0);
int row=i*stride;
int col=j*stride;
va=vld1q_f32(&input[row*in_size+col]);
vb=vld1q_f32(&input[(row+1)*in_size+col]);
vc=vld1q_f32(&input[(row+2)*in_size+col]);
va=vsetq_lane_f32(0,va,3);
vb=vsetq_lane_f32(0,vb,3);
vc=vsetq_lane_f32(0,vc,3);
max_v=vmaxq_f32(max_v,va);
max_v=vmaxq_f32(max_v,vb);
max_v=vmaxq_f32(max_v,vc);
float32x2_t temp_max_v=vmax_f32(vget_high_f32(max_v),vget_low_f32(max_v));
temp_max_v=vpmax_f32(temp_max_v,temp_max_v);
float max_val=vget_lane_f32(temp_max_v,0);
output[i*out_size+j]=max_val;
}
}
}
在不打开03优化时的运行时间为66ms：
在打开O3优化时的运行时间居然仅为6ms左右，简直震惊：
（3）cuda版本
使用cuda可以并行处理多个kernel的池化，每一个kernel对应的池化使用一个线程进行计算。具体代码如下：
__global__
void maxpool(float *input,const int in_size,const int ker_size,const int stride,float *output){
int row=blockIdx.y*blockDim.y+threadIdx.y;
int col=blockIdx.x*blockDim.x+threadIdx.x;
int out_size=(in_size-ker_size)/stride+1;
if(row<out_size&&col<out_size){
float max=0;
for(int i=0;i<ker_size;++i){
for(int j=0;j<ker_size;++j){
float curr=input[(row*stride+i)*in_size+(col*stride+j)];
max=max<curr?curr:max;
}
}
output[row*out_size+col]=max;
}
}
使用nvprof分析运行时间如下：
maxpool核函数运行时间为11.8ms，比neon实现居然慢了一倍，搞不懂，可能还能进行优化？
2.averagepool
（1）基础版本
基础版本和maxpool差不多，具体不在赘述。
void average_pool(const float *input,const int in_size,const int ker_size,const int stride,float *output){
int out_size=(in_size-ker_size)/stride+1;
for(int i=0;i<out_size;++i){
for(int j=0;j<out_size;++j){
for(int ker_i=0;ker_i<ker_size;++ker_i){
for(int ker_j=0;ker_j<ker_size;++ker_j){
int in_index=(i*stride+ker_i)*in_size+j*stride+ker_j;
int out_index=i*out_size+j;
output[out_index]=input[in_index];
}
}
}
}
}
不打开O3优化的具体运行结果如下：
打开O3优化的具体结果如下：
（2）neon版本
neon实现averagepool只需要将之前的maxpool算子中的求最大值操作改为求平均操作就可以了。具体参考代码如下：
void averagepool(const float *input,const int in_size,const int ker_size,const int stride,float *output){
int out_size=(in_size-ker_size)/stride+1;
for(int i=0;i<out_size;++i){
for(int j=0;j<out_size;++j){
int row=i*stride;
int col=j*stride;
float32x4_t va=vdupq_n_f32(0);
float32x4_t vb=vdupq_n_f32(0);
float32x4_t vc=vdupq_n_f32(0);
va=vld1q_f32(&input[(row*in_size)+col]);
vb=vld1q_f32(&input[(row+1)*in_size+col]);
vc=vld1q_f32(&input[(row+2)*in_size+col]);
va=vsetq_lane_f32(0,va,3);
vb=vsetq_lane_f32(0,vb,3);
vc=vsetq_lane_f32(0,vc,3);
float32x4_t tmp=vaddq_f32(vc,vaddq_f32(va,vb));
float32x2_t tmp_sum=vadd_f32(vget_low_f32(tmp),vget_high_f32(tmp));
float sum=vaddv_f32(tmp_sum);
float average=sum/(ker_size*ker_size);
output[i*out_size+j]=average;
}
}
}
不打开-O3是的运行时间为大约66ms：
打开-O3的运行时间为大约为9ms，再次感叹-O3的优化真的太给力了：
（3）cuda版本
使用cuda实现averagepool算子的步骤和最大池化是差不多的，但在每一个子矩阵的池化操作是将求取最大值改写为求取平均值，没有什么难度。具体代码如下：
__global__
void average_pool(const float *input,const int in_size,const int ker_size,const int stride,float *output){
int row=blockIdx.y*blockDim.y+threadIdx.y;
int col=blockIdx.x*blockDim.x+threadIdx.x;
int out_size=(in_size-ker_size)/stride+1;
if(row<out_size&&col<out_size){
float sum=0;
for(int ker_i=0;ker_i<ker_size;++ker_i){
for(int ker_j=0;ker_j<ker_size;++ker_j){
int index=(row*stride+ker_i)*in_size+col*stride+ker_j;
sum+=input[index];
}
}
float res=sum/(ker_size*ker_size);
output[row*out_size+col]=res;
}
}
针对640*640的输入矩阵的数据，具体执行时间为12ms左右：
三、im2col
1.cuda实现
不管是前面写过的卷积操作还是池化操作，其实都有一个类似的特点，也就是都将整个输入张量拆分成较小的子操作进行处理，比如卷积操作的处理尺度大多数为3x3、5x5等，池化操作也类似。但是在进行这两个操作的时候其实非常的耗时。从之前的两种算子在不同平台上（arm neon和cuda）的实现可以看出，某些算子在arm上利用neon指令处理得更快，原因其实是因为并没有充分利用到gpu上的计算核心，并且在进行卷积操作时，同一个子矩阵中的数据并不是连续的，因此无法充分利用好一个warp的并行，所以为了能够充分利用上gpu的并行能力，可以使用im2col的想法进行优化。
im2col的原理主要就是将输入矩阵中，需要进行操作的每一个子矩阵转变为一个列向量，然后使用矩阵乘法来进行计算，这样做得好处其实就是将一个子矩阵中的数据变为一个连续的列向量，更加符合cuda的处理逻辑。
下面实现一下im2col算子，平台是jetson nano，搭载一块maxwell架构gpu和arm架构cpu。使用测试用例为640*640的张量。3*3的kernel_size。
下面的写法是更能够充分利用gpu的计算核心，比上面的写法具有更小的粒度，每一个核心负责转移一个数据，而并不是一个子矩阵的数据，这样做就会更快了。
__global__
void im2col(const float *input,const int in_size,const int ker_size,const int stride,float *output){
int idy=blockIdx.y*blockDim.y+threadIdx.y;
int idx=blockIdx.x*blockDim.x+threadIdx.x;
int out_size=(in_size-ker_size)/stride+1;
int width=out_size*out_size;
int height=ker_size*ker_size;
if(idx<width && idy<height){
int row=idx/out_size;
int col=idx%out_size;
int input_row=row+(idy/ker_size);
int input_col=col+(idy%ker_size);
int input_index=input_row*in_size+input_col;
output[idy*width+idx]=input[input_index];
}
}
运行时间为3ms左右：
2.im2col加速卷积操作
在使用上面的im2col将输入矩阵按照需要卷积的子矩阵变换以后，就可以使用矩阵乘法进行并行计算了。
gpu上的每一个线程实现每一个子矩阵的乘加，具体代码如下：
void im2col(const float *input,const int in_size,const int ker_size,const int stride,float *output){
int idy=blockIdx.y*blockDim.y+threadIdx.y;
int idx=blockIdx.x*blockDim.x+threadIdx.x;
int out_size=(in_size-ker_size)/stride+1;
int width=out_size*out_size;
int height=ker_size*ker_size;
if(idx<width && idy<height){
int row=idx/out_size;
int col=idx%out_size;
int input_row=row+(idy/ker_size);
int input_col=col+(idy%ker_size);
int input_index=input_row*in_size+input_col;
output[idy*width+idx]=input[input_index];
}
}
__device__
void GEMM(const float *input,const int in_width,const int in_height,const float *kernel,const int ker_size,float *output){
int idx=blockIdx.y*blockDim.x+threadIdx.x;
if(idx<in_width){
float sum=0;
for(int i=0;i<in_height;++i){
sum+=input[i*in_width+idx]*kernel[i];
}
output[idx]=sum;
}
}
__global__
void conv2d(const float *input,const int in_size,const float *kernel,const int ker_size,const int stride,float *im2col_array,float *output){
const int out_size=(in_size-ker_size)/stride+1;
const int in_width=out_size*out_size;
const int in_height=ker_size*ker_size;
im2col(input,in_size,ker_size,stride,im2col_array);
GEMM(im2col_array,in_width,in_height,kernel,ker_size,output);
}
测试数据位640*640的矩阵，卷积核大小为3*3，步长为1。运行时间大致为9ms，比不适用im2col降低了4ms左右。
四、Reduce
可查看博客reduce
五、softmax
1.基础版本
softmax算子的原理其实就是根据各个元素在珍格格数组中的占比来计算出百分比，这个算子常常用在分类任务上。该算子的视线主要分为几个任务，首先是计算出整个数组中所有元素的最大值，其目的是用来将所有元素缩小到0-1之间，防止在指数计算的时候出现数据溢出。然后将所有元素计算其指数值，以及所有指数的和，最后一步是计算每一个元素的指数占指数和的比例。
具体代码实现如下：
void softmax(float *input,const int len,float *output){
float sum=0.0f;
float max_num=0.0f;
for(int i=0;i<len;++i){
max_num=max_num<input[i]?input[i]:max_num;
}
for(int i=0;i<len;i++){
input[i]=input[i]/max_num;
sum+=exp(input[i]);
}
for(int i=0;i<len;i++){
output[i]=exp(input[i])/sum;
}
}
使用100000个元素的数组进行测试，具体运行时间为22ms：
2.neon版本
在基础版本的基础上使用neon intrinsic函数进行加速，主要也是对三个阶段进行加速。首先是求最大值上的加速，这部分主要是使用neon寄存器4个元素4个元素的比较，找出最大值；第二部分是求每个元素的指数，并求所有元素的指数和，因为neon intrinsic没有exp函数直接进行指数运算，所以我使用的是指数函数的二级泰勒展开式进行指数的近似计算；第三步是计算每个元素在所有指数和的占比。另外还将循环展开，充分利用上-O3优化。
具体实现如下：
void softmax_neon(float *input,const int len,float *output){
float32x4_t max_vec=vdupq_n_f32(0); //求取input中的最大值
for(int i=0;i<len;i+=4){
float32x4_t tmp=vld1q_f32(&input[i]);
max_vec=vmaxq_f32(max_vec,tmp);
}
float32x2_t maxpair=vmax_f32(vget_low_f32(max_vec),vget_high_f32(max_vec));
maxpair=vmax_f32(maxpair,maxpair);
float max_val=vget_lane_f32(maxpair,0);
float32x4_t diver1=vdupq_n_f32(max_val); //所有元素除以最大值，防止数据溢出
float32x4_t sum_rig=vdupq_n_f32(0);
for(int i=0;i<len;i+=4){
float32x4_t tmp=vld1q_f32(&input[i]);
vst1q_f32(&input[i],vdivq_f32(tmp,diver1));
}
for(int i=0;i<len;i+=4){ //求解每一个元素的指数，并求出指数和
float32x4_t tmp=vld1q_f32(&input[i]);
float32x4_t one=vdupq_n_f32(1.0f);
float32x4_t a=vdupq_n_f32(0.99989f);
float32x4_t b=vdupq_n_f32(0.54030f);
float32x4_t res=one+tmp*(a+b*tmp);
sum_rig=vaddq_f32(sum_rig,res);
vst1q_f32(&input[i],res);
}
float32x2_t sumpair=vadd_f32(vget_low_f32(sum_rig),vget_high_f32(sum_rig));
sumpair=vadd_f32(sumpair,sumpair);
float sum_val=vget_lane_f32(sumpair,0);
float32x4_t diver2=vdupq_n_f32(sum_val);
for(int i=0;i<len;i+=4){
float32x4_t tmp=vld1q_f32(&input[i]);
float32x4_t result=vdivq_f32(tmp,diver2);
vst1q_f32(&output[i],result);
}
}
在使用指数函数的泰勒展开式求解指数时，对应的三个常量其实是根据经验分析所得的，可以根据任务的需求进行修改。另外在使用循环展开进行优化时，会使用到较多的寄存器，应该注意展开的程度。
使用neon加速运行时间为大概6ms左右：
打开-O3优化后的运行时间不到4ms：
3090
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
