第7章并行模式：卷积 - Stencil计算 (Parallel Patterns: Convolution - Stencil Computation)

一、导读

7.1 本章学习目标

通过本章的学习，读者将能够：

理解卷积（Convolution）和Stencil计算的基本概念及其在科学计算和信号处理中的广泛应用

掌握一维和二维卷积的并行化策略，能够设计高效的CUDA kernel

深入理解共享内存在卷积优化中的核心作用，学会使用halo/ghost区域减少全局内存访问

掌握边界处理的各种方法（零填充、镜像、周期性、常数），并能根据应用需求选择合适的方法

理解可分离卷积的原理，能够将高维卷积分解为多个一维卷积以提高效率

掌握常量内存和纹理内存在存储卷积核和优化访问模式中的应用

学会分析卷积kernel的性能瓶颈，包括计算强度、内存带宽和occupancy的权衡

了解高阶Stencil模式和变尺寸卷积核的处理方法

7.2 为什么本章如此重要

卷积是信号处理、图像处理、科学模拟和深度学习中最基础的运算之一。在图像处理中，模糊、锐化、边缘检测等操作本质上都是卷积；在科学计算中，有限差分法求解偏微分方程需要Stencil计算；在深度学习中，卷积神经网络（CNN）的核心就是卷积操作。

卷积的计算特点是每个输出元素依赖于输入数据的一个局部邻域，这种局部依赖性使得卷积天然适合并行化。然而，简单的并行实现往往效率低下，因为相邻线程会重复加载相同的全局内存数据。使用共享内存缓存复用数据可以将全局内存访问减少数倍，这是卷积优化的关键。

此外，卷积有多种变体：一维、二维、三维卷积；固定大小和可变大小卷积核；标准卷积和深度可分离卷积；全卷积和空洞卷积等。理解这些变体的特点和适用场景，对于在不同应用中选择合适的实现至关重要。

本章将从简单的一维卷积开始，逐步深入到复杂的二维和多维卷积，系统地介绍各种优化技术。通过详细的代码示例和性能分析，读者将掌握卷积并行化的核心技能，为后续学习更复杂的并行模式奠定基础。

7.3 核心问题引导

什么是卷积？它与相关（correlation）有什么区别？

一维卷积的并行化策略是什么？如何映射线程到输出元素？

为什么共享内存对卷积优化如此重要？halo区域的作用是什么？

如何处理卷积的边界条件？不同边界处理方法的应用场景是什么？

什么是可分离卷积？如何将二维卷积分解为两个一维卷积？

常量内存和纹理内存在卷积中各有什么优势？

卷积的计算强度（arithmetic intensity）是多少？它是计算受限还是内存受限？

大尺寸卷积核（如15×15或更大）应该如何优化？

二、核心概念详解

7.1 卷积与Stencil计算基础

卷积的定义：

卷积是一种数学运算，它将两个函数（通常是输入信号和滤波器/核）组合产生第三个函数。在离散情况下，一维卷积定义为：

$$(f * g)[n] = \sum_{k=-r}^{r} f[n-k] \cdot g[k+r]$$

其中，$f$是输入信号，$g$是卷积核（滤波器），$r$是核半径，核大小为$2r+1$。

在二维情况下，卷积定义为：

$$(F * G)[i,j] = \sum_{u=-r}^{r} \sum_{v=-r}^{r} F[i-u, j-v] \cdot G[u+r, v+r]$$

卷积 vs 相关：

卷积和相关（correlation）非常相似，区别在于卷积会将核翻转180度，而相关不会。在许多应用中（如图像处理），我们实际上执行的是相关操作，但习惯上仍称为"卷积"。在深度学习中，CNN的"卷积层"实际上执行的是相关运算。

Convolution:  output[i] = sum(input[i - k] * kernel[k])   // kernel flipped
Correlation:  output[i] = sum(input[i + k] * kernel[k])   // kernel not flipped

Stencil计算：

Stencil计算是卷积的推广，指每个输出点通过对输入点的某个固定模式（stencil）的邻域进行计算得到。Stencil可以是规则的（如矩形窗口）或不规则的（如自适应邻域）。有限差分法中的导数近似是典型的Stencil计算：

$$\frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x}$$

这是一个3点中心差分stencil。

应用领域：

图像处理：高斯模糊、Sobel边缘检测、Laplacian锐化等

科学模拟：热传导方程、波动方程、流体力学的有限差分离散

信号处理：低通滤波、高通滤波、特征提取

深度学习：卷积神经网络的卷积层

计算机视觉：特征检测、图像金字塔

7.2 一维卷积的并行化

基本思路：

一维卷积的并行化非常直观：每个输出元素可以由一个独立的线程计算。线程索引直接映射到输出数组的索引。

朴素实现：

cpp

// Naive 1D convolution using global memory only
__global__ void convolve1DNaive(const float* input, float* output, 
                                 const float* kernel, int N, int kernelRadius) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int kernelSize = 2 * kernelRadius + 1;
    
    if (idx < N) {
        float sum = 0.0f;
        
        for (int k = -kernelRadius; k <= kernelRadius; k++) {
            // Handle boundary: clamp to edge
            int neighborIdx = min(N - 1, max(0, idx + k));
            sum += input[neighborIdx] * kernel[k + kernelRadius];
        }
        
        output[idx] = sum;
    }
}

这个实现的缺点是每次卷积都需要从全局内存读取kernelSize个元素，相邻线程会重复读取相同的数据，导致大量的冗余全局内存访问。

性能分析：

对于一个大小为N的输入和大小为K的核：

全局内存读取次数：N × K（每个输出元素读取K个输入值）
全局内存写入次数：N
计算量：N × K次乘加

当K较大时，全局内存访问成为主要瓶颈。

7.3 使用共享内存优化一维卷积

优化思路：

利用共享内存缓存输入数据的局部窗口，使得每个输入元素只从全局内存读取一次，然后被多个线程复用。

Halo/Ghost区域：

为了计算Block边界处的输出元素，需要从相邻Block获取额外的输入数据。这些额外数据称为halo或ghost区域。对于一个半径为r的核，每个Block需要在两侧各扩展r个元素。

cpp

#define BLOCK_SIZE 256
#define KERNEL_RADIUS 4

__global__ void convolve1DShared(const float* input, float* output,
                                  const float* kernel, int N) {
    __shared__ float shared[BLOCK_SIZE + 2 * KERNEL_RADIUS];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * BLOCK_SIZE + tid;
    
    // Load main data into shared memory
    if (globalIdx < N) {
        shared[tid + KERNEL_RADIUS] = input[globalIdx];
    }
    
    // Load halo regions
    // Left halo
    if (tid < KERNEL_RADIUS) {
        int leftIdx = blockIdx.x * BLOCK_SIZE + tid - KERNEL_RADIUS;
        if (leftIdx >= 0) {
            shared[tid] = input[leftIdx];
        } else {
            shared[tid] = input[0];  // Boundary: clamp
        }
    }
    
    // Right halo
    if (tid < KERNEL_RADIUS) {
        int rightIdx = blockIdx.x * BLOCK_SIZE + BLOCK_SIZE + tid;
        if (rightIdx < N) {
            shared[BLOCK_SIZE + KERNEL_RADIUS + tid] = input[rightIdx];
        } else {
            shared[BLOCK_SIZE + KERNEL_RADIUS + tid] = input[N - 1];  // Clamp
        }
    }
    
    __syncthreads();
    
    // Compute convolution from shared memory
    if (globalIdx < N) {
        float sum = 0.0f;
        for (int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++) {
            sum += shared[tid + KERNEL_RADIUS + k] * kernel[k + KERNEL_RADIUS];
        }
        output[globalIdx] = sum;
    }
}

优化效果分析：

使用共享内存后：

全局内存读取次数：N + 2 × (N/BLOCK_SIZE) × KERNEL_RADIUS ≈ N（当BLOCK_SIZE远大于KERNEL_RADIUS时）
共享内存读取次数：N × K（但共享内存速度远快于全局内存）

全局内存访问减少了约K倍，这对于较大的核尤其显著。

资源使用：

共享内存：(BLOCK_SIZE + 2 * KERNEL_RADIUS) × sizeof(float)字节
对于BLOCK_SIZE=256，KERNEL_RADIUS=4：(256 + 8) × 4 = 1056字节

这远小于SM的共享内存容量，不会显著影响occupancy。

7.4 二维卷积的实现

二维卷积的挑战：

二维卷积比一维复杂，因为：

数据是二维的，索引计算更复杂

共享内存需要二维数组

Halo区域需要在两个维度上扩展

内存访问模式需要考虑合并访问

朴素二维卷积：

cpp

__global__ void convolve2DNaive(const float* input, float* output,
                                 const float* kernel, 
                                 int width, int height, int kernelRadius) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int kernelSize = 2 * kernelRadius + 1;
    
    if (col < width && row < height) {
        float sum = 0.0f;
        
        for (int kr = -kernelRadius; kr <= kernelRadius; kr++) {
            for (int kc = -kernelRadius; kc <= kernelRadius; kc++) {
                int r = min(height - 1, max(0, row + kr));
                int c = min(width - 1, max(0, col + kc));
                
                int kernelIdx = (kr + kernelRadius) * kernelSize + (kc + kernelRadius);
                sum += input[r * width + c] * kernel[kernelIdx];
            }
        }
        
        output[row * width + col] = sum;
    }
}

这个实现的问题是每个线程独立从全局内存读取kernelSize²个元素，大量重复访问。

使用共享内存的二维卷积：

cpp

#define TILE_SIZE 16
#define KERNEL_RADIUS 2

__global__ void convolve2DShared(const float* input, float* output,
                                  const float* kernel,
                                  int width, int height) {
    // Shared memory with halo region
    __shared__ float shared[TILE_SIZE + 2 * KERNEL_RADIUS]
                                 [TILE_SIZE + 2 * KERNEL_RADIUS];
    
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int baseCol = blockIdx.x * TILE_SIZE;
    int baseRow = blockIdx.y * TILE_SIZE;
    
    int col = baseCol + tx;
    int row = baseRow + ty;
    
    // Load input into shared memory with halo
    // Each thread may load multiple elements to cover the extended region
    
    // Main tile loading
    if (col < width && row < height) {
        shared[ty + KERNEL_RADIUS][tx + KERNEL_RADIUS] = input[row * width + col];
    }
    
    // Load halo regions
    // Top halo
    if (ty < KERNEL_RADIUS) {
        int topRow = baseRow + ty - KERNEL_RADIUS;
        if (topRow >= 0 && col < width) {
            shared[ty][tx + KERNEL_RADIUS] = input[topRow * width + col];
        } else if (col < width) {
            shared[ty][tx + KERNEL_RADIUS] = input[0 * width + col];  // Clamp
        }
    }
    
    // Bottom halo
    if (ty < KERNEL_RADIUS) {
        int bottomRow = baseRow + TILE_SIZE + ty;
        if (bottomRow < height && col < width) {
            shared[TILE_SIZE + KERNEL_RADIUS + ty][tx + KERNEL_RADIUS] = 
                input[bottomRow * width + col];
        } else if (col < width) {
            shared[TILE_SIZE + KERNEL_RADIUS + ty][tx + KERNEL_RADIUS] = 
                input[(height - 1) * width + col];  // Clamp
        }
    }
    
    // Left halo
    if (tx < KERNEL_RADIUS) {
        int leftCol = baseCol + tx - KERNEL_RADIUS;
        if (leftCol >= 0 && row < height) {
            shared[ty + KERNEL_RADIUS][tx] = input[row * width + leftCol];
        } else if (row < height) {
            shared[ty + KERNEL_RADIUS][tx] = input[row * width + 0];  // Clamp
        }
    }
    
    // Right halo
    if (tx < KERNEL_RADIUS) {
        int rightCol = baseCol + TILE_SIZE + tx;
        if (rightCol < width && row < height) {
            shared[ty + KERNEL_RADIUS][TILE_SIZE + KERNEL_RADIUS + tx] = 
                input[row * width + rightCol];
        } else if (row < height) {
            shared[ty + KERNEL_RADIUS][TILE_SIZE + KERNEL_RADIUS + tx] = 
                input[row * width + (width - 1)];  // Clamp
        }
    }
    
    // Corner halos (simplified: let boundary threads handle)
    __syncthreads();
    
    // Compute convolution from shared memory
    if (col < width && row < height) {
        float sum = 0.0f;
        int kernelSize = 2 * KERNEL_RADIUS + 1;
        
        for (int kr = 0; kr < kernelSize; kr++) {
            for (int kc = 0; kc < kernelSize; kc++) {
                sum += shared[ty + kr][tx + kc] * kernel[kr * kernelSize + kc];
            }
        }
        
        output[row * width + col] = sum;
    }
}

简化加载策略：

上面的代码分别加载各个halo区域，较为繁琐。一个更简洁的方法是让每个线程负责加载多个元素：

cpp

__global__ void convolve2DSharedSimple(const float* input, float* output,
                                        const float* kernel,
                                        int width, int height) {
    __shared__ float shared[TILE_SIZE + 2 * KERNEL_RADIUS]
                                 [TILE_SIZE + 2 * KERNEL_RADIUS];
    
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int baseCol = blockIdx.x * TILE_SIZE;
    int baseRow = blockIdx.y * TILE_SIZE;
    
    // Each thread loads multiple elements
    int extendedWidth = TILE_SIZE + 2 * KERNEL_RADIUS;
    int extendedHeight = TILE_SIZE + 2 * KERNEL_RADIUS;
    
    for (int ly = ty; ly < extendedHeight; ly += blockDim.y) {
        for (int lx = tx; lx < extendedWidth; lx += blockDim.x) {
            int globalCol = baseCol + lx - KERNEL_RADIUS;
            int globalRow = baseRow + ly - KERNEL_RADIUS;
            
            // Clamp boundary
            globalCol = min(width - 1, max(0, globalCol));
            globalRow = min(height - 1, max(0, globalRow));
            
            shared[ly][lx] = input[globalRow * width + globalCol];
        }
    }
    
    __syncthreads();
    
    // Compute convolution
    int col = baseCol + tx;
    int row = baseRow + ty;
    
    if (col < width && row < height) {
        float sum = 0.0f;
        int kernelSize = 2 * KERNEL_RADIUS + 1;
        
        for (int kr = 0; kr < kernelSize; kr++) {
            for (int kc = 0; kc < kernelSize; kc++) {
                sum += shared[ty + kr][tx + kc] * kernel[kr * kernelSize + kc];
            }
        }
        
        output[row * width + col] = sum;
    }
}

这种方法更简洁，且能正确处理所有halo区域，包括角落。

7.5 边界处理方法

卷积在边界处需要特殊处理，因为核的一部分可能超出输入范围。常见的边界处理方法有：

1. 零填充（Zero Padding）：

超出边界的位置视为0。这是最简单的方法，但会在边界处引入暗边。

cpp

float getValue(int r, int c, int height, int width) {
    if (r < 0 || r >= height || c < 0 || c >= width) {
        return 0.0f;
    }
    return input[r * width + c];
}

2. 边缘克隆（Clamp/Replicate）：

超出边界的位置使用最近的边界值。这避免了暗边问题，是最常用的方法。

cpp

int clampedRow = min(height - 1, max(0, r));
int clampedCol = min(width - 1, max(0, c));
return input[clampedRow * width + clampedCol];

3. 周期性边界（Wrap/Periodic）：

将输入视为周期性的，超出部分从另一侧绕回。适用于周期性信号或纹理。

cpp

int wrappedRow = ((r % height) + height) % height;
int wrappedCol = ((c % width) + width) % width;
return input[wrappedRow * width + wrappedCol];

4. 镜像边界（Mirror/Reflect）：

超出边界的部分以边界为轴镜像反射。这在某些信号处理应用中能减少边界伪影。

cpp

int mirroredRow = r;
if (mirroredRow < 0) mirroredRow = -mirroredRow;
if (mirroredRow >= height) mirroredRow = 2 * (height - 1) - mirroredRow;
mirroredRow = min(height - 1, max(0, mirroredRow));

5. 常数填充：

使用用户指定的常数值填充边界。

在GPU实现中，边界处理可以在加载共享内存时完成（如上例所示），也可以在计算时动态判断。前者效率更高，因为边界检查只做一次；后者更灵活，但每次访问都需要判断。

7.6 可分离卷积

原理：

某些二维核可以分解为两个一维核的外积，这样的核称为可分离核。例如，高斯核是可分离的：

$$G(x,y) = G_x(x) \cdot G_y(y)$$

如果核$K$可以分解为$K = u \cdot v^T$，则二维卷积可以分解为两次一维卷积：

$$I K = (I u) * v$$

其中第一次卷积沿行方向进行，第二次沿列方向进行。

性能优势：

对于$K \times K$的核：

标准二维卷积：每像素$K^2$次乘加
可分离卷积：每像素$2K$次乘加

对于$K=5$，标准方法需要25次乘加，可分离方法只需10次，减少60%。

实现：

cpp

// First pass: horizontal convolution
__global__ void convolveHorizontal(const float* input, float* tempOutput,
                                    const float* kernelX,
                                    int width, int height, int kernelRadius) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col < width && row < height) {
        float sum = 0.0f;
        
        for (int k = -kernelRadius; k <= kernelRadius; k++) {
            int c = min(width - 1, max(0, col + k));
            sum += input[row * width + c] * kernelX[k + kernelRadius];
        }
        
        tempOutput[row * width + col] = sum;
    }
}

// Second pass: vertical convolution
__global__ void convolveVertical(const float* tempInput, float* output,
                                  const float* kernelY,
                                  int width, int height, int kernelRadius) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col < width && row < height) {
        float sum = 0.0f;
        
        for (int k = -kernelRadius; k <= kernelRadius; k++) {
            int r = min(height - 1, max(0, row + k));
            sum += tempInput[r * width + col] * kernelY[k + kernelRadius];
        }
        
        output[row * width + col] = sum;
    }
}

中间存储：

可分离卷积需要一个临时缓冲区存储第一次卷积的结果。对于大图像，这可能占用显著的显存。可以使用CUDA流重叠两次卷积的执行，或者使用双缓冲技术。

判断核是否可分离：

不是所有核都可分离。判断方法是对核矩阵进行奇异值分解（SVD），如果只有一个非零奇异值，则核可分离。在实际应用中，可以通过检查核的秩来判断。

7.7 常量内存与纹理内存优化

常量内存存储核：

卷积核通常是只读的常量数据，适合存储在常量内存中：

cpp

__constant__ float constKernel[25];  // For 5x5 kernel

// Copy kernel to constant memory
cudaMemcpyToSymbol(constKernel, h_kernel, kernelSize * sizeof(float));

// Use in kernel
__global__ void convolveWithConstantKernel(...) {
    // Access constKernel[i] - benefits from broadcast when all threads read same index
}

当warp内所有线程读取相同的核元素时，常量内存的广播机制使得访问极其高效。

纹理内存用于输入：

对于具有空间局部性的卷积，纹理内存可以提供额外的缓存优势：

cpp

texture<float, 2, cudaReadModeElementType> texInput;

// Bind input to texture
cudaBindTextureToArray(texInput, cuArray, channelDesc);

// Fetch in kernel
float val = tex2D(texInput, col, row);

纹理缓存针对二维空间局部性优化，对于卷积这种规则的空间访问模式特别有效。

7.8 性能分析与优化策略

计算强度：

计算强度（Arithmetic Intensity）定义为每字节内存传输的浮点运算次数。对于卷积：

计算量：$N \times K$次乘加（一维）或$N \times K^2$（二维）
内存传输：读取输入$N$个元素，写入输出$N$个元素，读取核$K$个元素

对于大$K$，计算强度约为$K/2$（一维）或$K^2/2$（二维）。当$K$较小时，卷积通常是内存带宽受限；当$K$较大时，逐渐变为计算受限。

优化策略总结：

小核（K ≤ 7）：使用共享内存，确保合并访问，核放入常量内存

中等核（7 < K ≤ 15）：使用共享内存，考虑可分离分解

大核（K > 15）：必须使用可分离分解（如果可能），或使用FFT-based卷积

多通道：对于RGB等多通道数据，可以考虑通道间的并行化

批量处理：对于多幅图像的卷积，使用批处理提高GPU利用率

三、程序员的视角

3.1 在实际项目中的应用

图像滤镜：Photoshop等软件中的模糊、锐化、浮雕等效果都是卷积操作。实时预览需要高效的GPU实现。

边缘检测：Canny、Sobel、Prewitt等边缘检测算子本质上是特定核的卷积。在自动驾驶和机器人视觉中，边缘检测是预处理的关键步骤。

深度学习：CNN的卷积层是现代计算机视觉的基础。虽然cuDNN提供了高度优化的实现，但理解底层原理有助于调试和优化网络。

科学模拟：热传导、扩散方程的数值求解使用Stencil计算。在气候模拟和材料科学中，这些计算占据大量时间。

3.2 常见陷阱与误区

忘记同步：在加载共享内存后忘记调用__syncthreads()，导致部分线程使用未初始化的数据。

Bank Conflict：二维共享内存访问可能导致bank conflict。使用填充（如[TILE_SIZE+1][TILE_SIZE+1]）可以避免。

边界处理错误：没有正确处理边界，导致访问越界或结果错误。始终在加载时进行边界检查。

核不可分离却强行分解：不是所有核都可分离。分解前需要验证核的可分离性。

过度优化小核：对于3×3这样的小核，共享内存优化的收益有限，因为数据复用率不高。有时直接使用全局内存配合缓存就足够了。

3.3 与其他层级的联系

与内存层次的联系：卷积是展示内存层次优化价值的经典案例。从全局内存到共享内存再到常量内存，每一层都发挥其优势。

与算法设计的联系：可分离卷积展示了算法层面的优化如何带来数量级的性能提升。程序员应该始终寻找问题的结构特性。

与硬件的联系：不同GPU架构的共享内存大小和bank数量不同，最优的tile大小可能需要调整。

四、代码实战

4.1 完整的一维卷积对比

cpp

#include <cuda_runtime.h>
#include <stdio.h>

#define N 1000000
#define KERNEL_RADIUS 5
#define BLOCK_SIZE 256

// Naive version
__global__ void convolve1DNaive(const float* input, float* output,
                                 const float* kernel, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        float sum = 0.0f;
        for (int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++) {
            int neighborIdx = min(N - 1, max(0, idx + k));
            sum += input[neighborIdx] * kernel[k + KERNEL_RADIUS];
        }
        output[idx] = sum;
    }
}

// Shared memory version
__global__ void convolve1DShared(const float* input, float* output,
                                  const float* kernel, int N) {
    __shared__ float shared[BLOCK_SIZE + 2 * KERNEL_RADIUS];
    
    int tid = threadIdx.x;
    int globalIdx = blockIdx.x * BLOCK_SIZE + tid;
    
    // Load main data
    if (globalIdx < N) {
        shared[tid + KERNEL_RADIUS] = input[globalIdx];
    }
    
    // Load halos
    if (tid < KERNEL_RADIUS) {
        int leftIdx = blockIdx.x * BLOCK_SIZE + tid - KERNEL_RADIUS;
        shared[tid] = (leftIdx >= 0) ? input[leftIdx] : input[0];
        
        int rightIdx = blockIdx.x * BLOCK_SIZE + BLOCK_SIZE + tid;
        shared[BLOCK_SIZE + KERNEL_RADIUS + tid] = 
            (rightIdx < N) ? input[rightIdx] : input[N - 1];
    }
    
    __syncthreads();
    
    if (globalIdx < N) {
        float sum = 0.0f;
        for (int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++) {
            sum += shared[tid + KERNEL_RADIUS + k] * kernel[k + KERNEL_RADIUS];
        }
        output[globalIdx] = sum;
    }
}

void benchmark1DConvolution() {
    printf("=== 1D Convolution Benchmark ===\n");
    
    size_t dataSize = N * sizeof(float);
    size_t kernelSize = (2 * KERNEL_RADIUS + 1) * sizeof(float);
    
    float *h_input = (float*)malloc(dataSize);
    float *h_kernel = (float*)malloc(kernelSize);
    
    for (int i = 0; i < N; i++) h_input[i] = (float)i / N;
    for (int i = 0; i < 2 * KERNEL_RADIUS + 1; i++) h_kernel[i] = 1.0f / (2 * KERNEL_RADIUS + 1);
    
    float *d_input, *d_output, *d_kernel;
    cudaMalloc((void**)&d_input, dataSize);
    cudaMalloc((void**)&d_output, dataSize);
    cudaMalloc((void**)&d_kernel, kernelSize);
    
    cudaMemcpy(d_input, h_input, dataSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_kernel, h_kernel, kernelSize, cudaMemcpyHostToDevice);
    
    int gridSize = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // Naive
    cudaEventRecord(start);
    for (int i = 0; i < 50; i++) {
        convolve1DNaive<<<gridSize, BLOCK_SIZE>>>(d_input, d_output, d_kernel, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float msNaive = 0;
    cudaEventElapsedTime(&msNaive, start, stop);
    
    // Shared
    cudaEventRecord(start);
    for (int i = 0; i < 50; i++) {
        convolve1DShared<<<gridSize, BLOCK_SIZE>>>(d_input, d_output, d_kernel, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float msShared = 0;
    cudaEventElapsedTime(&msShared, start, stop);
    
    printf("Naive:     %.3f ms\n", msNaive / 50);
    printf("Shared:    %.3f ms\n", msShared / 50);
    printf("Speedup:   %.2fx\n", msNaive / msShared);
    
    cudaFree(d_input);
    cudaFree(d_output);
    cudaFree(d_kernel);
    free(h_input);
    free(h_kernel);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

int main() {
    benchmark1DConvolution();
    return 0;
}

4.2 二维可分离高斯卷积

cpp

#include <cuda_runtime.h>
#include <stdio.h>
#include <cmath>

#define WIDTH 2048
#define HEIGHT 2048
#define GAUSSIAN_RADIUS 3
#define TILE_SIZE 16

// Create 1D Gaussian kernel
void createGaussianKernel(float* kernel, int radius, float sigma) {
    float sum = 0.0f;
    for (int i = -radius; i <= radius; i++) {
        kernel[i + radius] = expf(-(i * i) / (2 * sigma * sigma));
        sum += kernel[i + radius];
    }
    for (int i = 0; i < 2 * radius + 1; i++) {
        kernel[i] /= sum;
    }
}

// Horizontal pass
__global__ void gaussianHorizontal(const float* input, float* output,
                                    const float* kernel,
                                    int width, int height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col < width && row < height) {
        float sum = 0.0f;
        for (int k = -GAUSSIAN_RADIUS; k <= GAUSSIAN_RADIUS; k++) {
            int c = min(width - 1, max(0, col + k));
            sum += input[row * width + c] * kernel[k + GAUSSIAN_RADIUS];
        }
        output[row * width + col] = sum;
    }
}

// Vertical pass
__global__ void gaussianVertical(const float* input, float* output,
                                  const float* kernel,
                                  int width, int height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col < width && row < height) {
        float sum = 0.0f;
        for (int k = -GAUSSIAN_RADIUS; k <= GAUSSIAN_RADIUS; k++) {
            int r = min(height - 1, max(0, row + k));
            sum += input[r * width + col] * kernel[k + GAUSSIAN_RADIUS];
        }
        output[row * width + col] = sum;
    }
}

void benchmarkSeparableGaussian() {
    printf("=== Separable Gaussian Convolution (%dx%d) ===\n", WIDTH, HEIGHT);
    
    size_t imageSize = WIDTH * HEIGHT * sizeof(float);
    size_t kernelSize = (2 * GAUSSIAN_RADIUS + 1) * sizeof(float);
    
    float *h_input = (float*)malloc(imageSize);
    for (int i = 0; i < WIDTH * HEIGHT; i++) {
        h_input[i] = (float)(i % 256);
    }
    
    float h_kernel[2 * GAUSSIAN_RADIUS + 1];
    createGaussianKernel(h_kernel, GAUSSIAN_RADIUS, 1.0f);
    
    float *d_input, *d_temp, *d_output, *d_kernel;
    cudaMalloc((void**)&d_input, imageSize);
    cudaMalloc((void**)&d_temp, imageSize);
    cudaMalloc((void**)&d_output, imageSize);
    cudaMalloc((void**)&d_kernel, kernelSize);
    
    cudaMemcpy(d_input, h_input, imageSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_kernel, h_kernel, kernelSize, cudaMemcpyHostToDevice);
    
    dim3 blockSize(TILE_SIZE, TILE_SIZE);
    dim3 gridSize((WIDTH + TILE_SIZE - 1) / TILE_SIZE, 
                  (HEIGHT + TILE_SIZE - 1) / TILE_SIZE);
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    cudaEventRecord(start);
    for (int i = 0; i < 20; i++) {
        gaussianHorizontal<<<gridSize, blockSize>>>(d_input, d_temp, d_kernel, 
                                                     WIDTH, HEIGHT);
        gaussianVertical<<<gridSize, blockSize>>>(d_temp, d_output, d_kernel,
                                                   WIDTH, HEIGHT);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms = 0;
    cudaEventElapsedTime(&ms, start, stop);
    
    printf("Separable Gaussian: %.3f ms per iteration\n", ms / 20);
    
    cudaFree(d_input);
    cudaFree(d_temp);
    cudaFree(d_output);
    cudaFree(d_kernel);
    free(h_input);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

int main() {
    benchmarkSeparableGaussian();
    return 0;
}

五、深入思考

5.1 历史背景

卷积的概念最早由数学家D'Alembert在18世纪提出，后来在信号处理和图像处理中得到广泛应用。在计算机图形学中，卷积是图像滤波的基础操作。

GPU加速卷积的历史可以追溯到早期的GPGPU研究，当时研究者通过将图像渲染为纹理来间接执行卷积。CUDA的出现使得直接编程卷积成为可能。cuDNN库的发布进一步推动了卷积在深度学习中的应用，其高度优化的实现达到了接近硬件理论峰值的性能。

近年来，Winograd算法和FFT-based卷积为大核卷积提供了新的优化方向。Tensor Core的引入使得混合精度卷积成为可能，进一步提升了吞吐量。

5.2 前沿方向

Winograd卷积：通过数学变换减少乘法次数，特别适合小核（3×3, 5×5）。Winograd已被cuDNN采用。

FFT卷积：对于大核，在频域执行卷积更高效。复杂度从O(K²)降为O(log K)。

稀疏卷积：利用输入或核的稀疏性跳过零值计算。在点云处理和稀疏体素渲染中有应用。

动态卷积：核权重根据输入动态生成，用于注意力机制和自适应滤波。

5.3 延伸阅读

"Real-Time Convolution on GPUs", NVIDIA Developer Blog

cuDNN Documentation - Convolution Operations

Lavin et al., "Fast Algorithms for Convolutional Neural Networks", CVPR 2016

MathWorks, "Implementing Convolution on GPU"

六、实验笔记

实验1：共享内存优化效果量化

目的：测量共享内存对卷积性能的改善。

步骤：

实现朴素和共享内存版本的一维卷积

测试不同核大小（3, 5, 7, 9, 11）下的性能

绘制加速比随核大小的变化曲线

预期：核越大，共享内存的加速效果越明显，因为数据复用率更高。

实验2：可分离vs标准二维卷积

目的：验证可分离卷积的性能优势。

步骤：

实现标准和可分离版本的5×5高斯卷积

比较执行时间和数值误差

分析两种方法的内存访问模式

预期：可分离版本应快约2倍，数值误差在可接受范围内。

七、本章小结

本章详细介绍了卷积这一重要的并行模式。我们从一维卷积入手，展示了如何使用共享内存和halo区域优化全局内存访问。然后扩展到二维卷积，讨论了更复杂的索引计算和边界处理。

可分离卷积是一个重要的优化技术，它将高维卷积分解为多个一维卷积，显著减少了计算量。但需要注意，并非所有核都可分离，使用前需要验证。

边界处理是卷积实现中容易被忽视但至关重要的部分。不同的边界处理方法适用于不同的应用场景，程序员需要根据具体需求选择。

常量内存和纹理内存为卷积提供了额外的优化手段。常量内存的广播机制适合存储核，纹理缓存的空间局部性优化适合输入数据。

通过代码示例，我们实践了一维和二维卷积的实现，以及可分离高斯卷积的优化。这些技能可以直接应用于图像处理、科学计算和深度学习等领域。

八、习题

卷积定义：解释卷积与相关的区别。在什么情况下两者等价？

一维索引计算：对于长度为N的输入和半径为r的核，输出长度是多少？如果使用"valid"卷积（不处理边界），输出长度又是多少？

共享内存大小：对于TILE_SIZE=16和KERNEL_RADIUS=3的二维卷积，共享内存需要多少字节？如果改为TILE_SIZE=32，共享内存需求如何变化？

可分离性判断：以下3×3核是否可分离？

```

[1 2 1]

[2 4 2]