我們之前介紹了cudnn調用api直接實現卷積，本文我們探究手動實現。
對于直接使用for循環在cpu上的實現方法，就不過多介紹，只要了解卷積的原理，就很容易實現。
im2col 的核心思想
im2col = image to column
把輸入 feature map 的每個卷積感受野（sliding window）展開成一列向量
卷積核也展開成一個行向量
然后把卷積轉化成 矩陣乘法（GEMM, General Matrix Multiply）
舉例
假設：
輸入：1 channel, 4×4
卷積核：1×2×2
步長 stride=1
輸入：
1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16
卷積核 2×2 的每個 sliding window：
[[1,2],[5,6]] → 展開為 [1,2,5,6]
[[2,3],[6,7]] → 展開為 [2,3,6,7]
...
用 im2col 后：
X_col = [[1, 2, 5, 6], # 第一個位置[2, 3, 6, 7], # 第二個位置[3, 4, 7, 8],...]
卷積核也展開成：
W_col = [w1, w2, w3, w4]
然后 卷積計算就變成矩陣乘法：
輸出每個位置就是矩陣乘法的一個元素
對多通道、多卷積核也可以批量做
說人話，就是把卷積核每次對準的這塊矩陣區域展平成向量；比如X_col的第一行，與W_col作向量乘法，結果就是第一次卷積得到的結果。
但是一般我們會將W_col的參數作為一行，所以X_col實際存儲需要轉置一下，這樣才符合矩陣乘法的要求
為什么效率高？
GEMM 有高度優化
BLAS/cuBLAS/cuDNN 都對矩陣乘法做了很多優化
可以充分利用 SIMD / GPU 并行
循環嵌套少
原始卷積是 6 層循環（batch, output channel, height, width, input channel, kernel height/width）
im2col + GEMM → 只要一次矩陣乘法
容易擴展到多通道、多 batch
代碼實現
代碼中對應的X_col就是按列存儲展開的元素
#ifndef __CUDACC__
#define __CUDACC__
#endif
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cudnn.h>
#include <cublas_v2.h>
#include <iostream>
#include<cstdio>
#include <cmath>
#include <cstdlib>
#include <vector>__global__ void im2col_kernel(const float *data_im,//輸入圖片數據 [C,H,W]int channels,int height,int width,// 輸入通道數 C,輸入高度 H, 寬度 Wint ksize,int pad,int stride,// 卷積核大小 (假設方形: kH=kW=ksize),padding 大小,卷積步長int height_col,int width_col,// 輸出特征圖的高寬float *data_col){ // im2col 展開的矩陣輸出 [C*ksize*ksize, H_out*W_out]int index=blockIdx.x*blockDim.x+threadIdx.x;//index: 每個線程負責展開一個卷積窗口中的一個元素。int n=channels*ksize*ksize*height_col*width_col;//n: 總的元素數if(index<n){//對比一下w_out*h_out*k_idx*k_w*k_h與n就知道為什么要這么計算這五個變量了//w_out, h_out: 代表當前處理的輸出特征圖位置 (卷積結果的坐標)。int w_out=index%width_col;int h_out=(index/width_col)%height_col;//k_idx → (c_in, k_h, k_w): 把卷積核的 index 分解成通道號、核內行列。int k_idx=(index/width_col/height_col);int k_w=k_idx%ksize;int k_h=(k_idx/ksize)%ksize;//c_in輸入通道索引int c_in=k_idx/(ksize*ksize);//im_row, im_col: 對應輸入圖片上的真實位置（考慮 stride、padding 后）//這是經典的把輸出坐標映射回輸入坐標的公式int im_row=h_out*stride-pad+k_h;int im_col=w_out*stride-pad+k_w;//col_index: data_col 的存儲索引int col_index=(c_in*ksize*ksize+k_h*ksize+k_w)* (height_col * width_col) + h_out * width_col + w_out;float val=0;//如果 im_row, im_col 在輸入圖像范圍內 → 取值；否則屬于padding → 0if(im_row>=0&&im_row<height&&im_col>=0&&im_col<width){val=data_im[(c_in*height+im_row)*width+im_col];}data_col[col_index]=val;}
}
void im2col_gpu(const float *d_im,int channels,int height,int width,int ksize,int pad,int stride,float *d_col){//計算輸出大小：就是卷積輸出 H_out, W_out 的公式。int height_col=(height +2*pad-ksize)/stride+1;int width_col=(width+2*pad-ksize)/stride+1;int n=channels*ksize*ksize*height_col*width_col;int threads=256;int blocks=(n+threads-1)/threads;im2col_kernel<<<blocks,threads>>>(d_im,channels,height,width,ksize,pad,stride,height_col,width_col,d_col);cudaDeviceSynchronize();
}
void conv_forward_im2col(const float *d_input,// 輸入圖像 [C,H,W]const float *d_weight,// 卷積核 [K,C,ksize,ksize]float *d_output,//輸出 [K,H_out,W_out]int C,int H,int W,// 輸入通道數，高，寬int K,// 卷積核數量 (輸出通道數)int ksize,int stride,int pad,// 核大小，步長，填充cublasHandle_t &handle){int H_out=(H+2*pad-ksize)/stride+1;int W_out=(W+2*pad-ksize)/stride+1;float *d_col;// im2col buffer: [C*ksize*ksize, H_out*W_out]cudaMalloc(&d_col,sizeof(float)*C*ksize*ksize*H_out*W_out);im2col_gpu(d_input,C,H,W,ksize,pad,stride,d_col);//原本是W*x的順序，但是由于cublasSgemm函數的特性，寫的時候參照實際傳遞的方式// GEMM: d_weight:[K, C*ksize*ksize] * d_col:[C*ksize*ksize, H_out*W_out] = [K, H_out*W_out]const float alpha=1.0f,beta=0.0f;cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,H_out*W_out,K,C*ksize*ksize,&alpha,d_col,H_out*W_out,d_weight,C*ksize*ksize,&beta,d_output,H_out*W_out);cudaFree(d_col);
}
int main() {int C=1, H=5, W=5;int K=1, ksize=3, stride=1, pad=1;std::vector<float> h_input(C*H*W, 1.0f); // 全1輸入std::vector<float> h_weight(K*C*ksize*ksize, 1.0f); // 全1卷積核std::vector<float> h_output(K*H*W);float *d_input, *d_weight, *d_output;cudaMalloc(&d_input, h_input.size()*sizeof(float));cudaMalloc(&d_weight, h_weight.size()*sizeof(float));cudaMalloc(&d_output, h_output.size()*sizeof(float));cudaMemcpy(d_input, h_input.data(), h_input.size()*sizeof(float), cudaMemcpyHostToDevice);cudaMemcpy(d_weight, h_weight.data(), h_weight.size()*sizeof(float), cudaMemcpyHostToDevice);cublasHandle_t handle;cublasCreate(&handle);conv_forward_im2col(d_input, d_weight, d_output,C,H,W,K,ksize,stride,pad, handle);cudaMemcpy(h_output.data(), d_output, h_output.size()*sizeof(float), cudaMemcpyDeviceToHost);cublasDestroy(handle);// 打印結果int H_out=(H+2*pad-ksize)/stride+1;int W_out=(W+2*pad-ksize)/stride+1;std::cout << "Output (" << K << "," << H_out << "," << W_out << "):\n";for(int i=0;i<K;i++){for(int h=0;h<H_out;h++){for(int w=0;w<W_out;w++){std::cout << h_output[i*H_out*W_out+h*W_out+w] << " ";}std::cout << "\n";}}cudaFree(d_input);cudaFree(d_weight);cudaFree(d_output);return 0;
}
