以下是对“stencil-convolution”技能的中文翻译：

---

name: stencil-convolution description: 针对GPU优化的模板和卷积模式实现的专家技能。设计带有光环的平铺模板算法，实现2D/3D卷积核，优化边界条件处理，应用时间阻塞技术，生成可分离的滤波器实现，并分析模板内存带宽。 allowed-tools: Bash(*) 读写编辑Glob Grep WebFetch metadata: author: babysitter-sdk version: “1.0.0” category: domain-algorithms backlog-id: SK-013

stencil-convolution

你是 stencil-convolution - 一个针对GPU优化的模板和卷积模式实现的专家技能。这项技能为科学计算、图像处理和需要邻域计算的数值模拟提供了AI驱动的模板和卷积操作，包括：

• 设计带有光环的平铺模板算法
• 实现2D/3D卷积核
• 优化边界条件处理
• 应用时间阻塞技术
• 生成可分离的滤波器实现
• 配置共享内存平铺策略
• 分析模板内存带宽
• 支持多分辨率模板

概述

这项技能使得AI驱动的模板和卷积操作成为可能，包括：

• 设计带有光环的平铺模板算法
• 实现2D/3D卷积核
• 优化边界条件处理
• 应用时间阻塞技术
• 生成可分离的滤波器实现
• 配置共享内存平铺策略
• 分析模板内存带宽
• 支持多分辨率模板

先决条件

• NVIDIA CUDA Toolkit 11.0+
• 计算能力3.5+的GPU
• 理解内存合并模式
• Nsight Compute用于内存分析
• 可选：cuDNN用于优化卷积

能力

1. 基本2D模板（5点拉普拉斯）

// 简单的5点模板（用于比较）
__global__ void laplacian2D_naive(
    float* out, const float* in,
    int width, int height
) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;

    if (x >= 1 && x < width - 1 && y >= 1 && y < height - 1) {
        int idx = y * width + x;
        out[idx] = -4.0f * in[idx]
                 + in[idx - 1]      // 左
                 + in[idx + 1]      // 右
                 + in[idx - width]  // 上
                 + in[idx + width]; // 下
    }
}

2. 带有共享内存和光环的平铺模板

#define TILE_X 32
#define TILE_Y 32
#define HALO 1

__global__ void laplacian2D_tiled(
    float* out, const float* in,
    int width, int height
) {
    // 带有光环的共享内存
    __shared__ float tile[TILE_Y + 2 * HALO][TILE_X + 2 * HALO];

    // 全局坐标
    int gx = blockIdx.x * TILE_X + threadIdx.x;
    int gy = blockIdx.y * TILE_Y + threadIdx.y;

    // 共享内存中的局部坐标（偏移光环）
    int lx = threadIdx.x + HALO;
    int ly = threadIdx.y + HALO;

    // 加载中心平铺
    if (gx < width && gy < height) {
        tile[ly][lx] = in[gy * width + gx];
    }

    // 加载光环区域
    // 左光环
    if (threadIdx.x < HALO && gx >= HALO) {
        tile[ly][lx - HALO] = in[gy * width + (gx - HALO)];
    }
    // 右光环
    if (threadIdx.x >= TILE_X - HALO && gx + HALO < width) {
        tile[ly][lx + HALO] = in[gy * width + (gx + HALO)];
    }
    // 顶部光环
    if (threadIdx.y < HALO && gy >= HALO) {
        tile[ly - HALO][lx] = in[(gy - HALO) * width + gx];
    }
    // 底部光环
    if (threadIdx.y >= TILE_Y - HALO && gy + HALO < height) {
        tile[ly + HALO][lx] = in[(gy + HALO) * width + gx];
    }

    // 角落光环（如果需要更大的模板）
    // ...

    __syncthreads();

    // 使用共享内存计算模板
    if (gx >= 1 && gx < width - 1 && gy >= 1 && gy < height - 1) {
        out[gy * width + gx] = -4.0f * tile[ly][lx]
                             + tile[ly][lx - 1]
                             + tile[ly][lx + 1]
                             + tile[ly - 1][lx]
                             + tile[ly + 1][lx];
    }
}

3. 可配置半径的通用N点模板

template <int RADIUS>
__global__ void stencil2D_generic(
    float* out, const float* in,
    const float* weights,  // 模板权重
    int width, int height
) {
    extern __shared__ float tile[];

    const int TILE_X = blockDim.x;
    const int TILE_Y = blockDim.y;
    const int TILE_PITCH = TILE_X + 2 * RADIUS;

    int gx = blockIdx.x * TILE_X + threadIdx.x;
    int gy = blockIdx.y * TILE_Y + threadIdx.y;

    int lx = threadIdx.x + RADIUS;
    int ly = threadIdx.y + RADIUS;

    // 加载平铺带光环
    // ...（类似于上面但通用）

    __syncthreads();

    if (gx >= RADIUS && gx < width - RADIUS &&
        gy >= RADIUS && gy < height - RADIUS) {

        float result = 0.0f;
        int wIdx = 0;

        // 应用模板权重
        for (int dy = -RADIUS; dy <= RADIUS; dy++) {
            for (int dx = -RADIUS; dx <= RADIUS; dx++) {
                result += weights[wIdx++] * tile[(ly + dy) * TILE_PITCH + (lx + dx)];
            }
        }

        out[gy * width + gx] = result;
    }
}

4. 任意核的2D卷积

#define CONV_TILE_X 32
#define CONV_TILE_Y 32
#define MAX_KERNEL_RADIUS 8

// 常量内存中的核权重，快速访问
__constant__ float c_kernel[(2 * MAX_KERNEL_RADIUS + 1) * (2 * MAX_KERNEL_RADIUS + 1)];

__global__ void convolution2D(
    float* out, const float* in,
    int width, int height,
    int kernelRadius
) {
    extern __shared__ float tile[];

    int TILE_PITCH = CONV_TILE_X + 2 * kernelRadius;

    int gx = blockIdx.x * CONV_TILE_X + threadIdx.x;
    int gy = blockIdx.y * CONV_TILE_Y + threadIdx.y;

    int lx = threadIdx.x + kernelRadius;
    int ly = threadIdx.y + kernelRadius;

    // 加载中心
    if (gx < width && gy < height) {
        tile[ly * TILE_PITCH + lx] = in[gy * width + gx];
    } else {
        tile[ly * TILE_PITCH + lx] = 0.0f;  // 零填充
    }

    // 加载光环与边界处理
    // 左
    if (threadIdx.x < kernelRadius) {
        int srcX = gx - kernelRadius;
        tile[ly * TILE_PITCH + (lx - kernelRadius)] =
            (srcX >= 0 && gy < height) ? in[gy * width + srcX] : 0.0f;
    }
    // 右
    if (threadIdx.x >= CONV_TILE_X - kernelRadius) {
        int srcX = gx + kernelRadius;
        tile[ly * TILE_PITCH + (lx + kernelRadius)] =
            (srcX < width && gy < height) ? in[gy * width + srcX] : 0.0f;
    }
    // 顶部和底部（类似模式）
    // ...

    __syncthreads();

    if (gx < width && gy < height) {
        float sum = 0.0f;
        int kernelSize = 2 * kernelRadius + 1;

        for (int ky = -kernelRadius; ky <= kernelRadius; ky++) {
            for (int kx = -kernelRadius; kx <= kernelRadius; kx++) {
                int kidx = (ky + kernelRadius) * kernelSize + (kx + kernelRadius);
                sum += c_kernel[kidx] * tile[(ly + ky) * TILE_PITCH + (lx + kx)];
            }
        }

        out[gy * width + gx] = sum;
    }
}

5. 可分离卷积（2次通过以提高性能）

// 可分离卷积更快：O(2*r) vs O(r^2)
// 第一次通过：水平卷积
__global__ void convolutionRow(
    float* out, const float* in,
    int width, int height, int radius
) {
    extern __shared__ float tile[];

    int gx = blockIdx.x * blockDim.x + threadIdx.x;
    int gy = blockIdx.y;

    int TILE_WIDTH = blockDim.x + 2 * radius;
    int lx = threadIdx.x + radius;

    // 加载带光环
    if (gx < width) {
        tile[lx] = in[gy * width + gx];
    }
    if (threadIdx.x < radius) {
        tile[threadIdx.x] = (gx >= radius) ? in[gy * width + gx - radius] : 0.0f;
        tile[lx + blockDim.x] = (gx + blockDim.x < width) ?
            in[gy * width + gx + blockDim.x] : 0.0f;
    }

    __syncthreads();

    if (gx < width) {
        float sum = 0.0f;
        for (int k = -radius; k <= radius; k++) {
            sum += c_kernelRow[k + radius] * tile[lx + k];
        }
        out[gy * width + gx] = sum;
    }
}

// 第二次通过：垂直卷积
__global__ void convolutionColumn(
    float* out, const float* in,
    int width, int height, int radius
) {
    extern __shared__ float tile[];

    int gx = blockIdx.x;
    int gy = blockIdx.y * blockDim.y + threadIdx.y;

    int TILE_HEIGHT = blockDim.y + 2 * radius;
    int ly = threadIdx.y + radius;

    // 加载带光环
    if (gy < height) {
        tile[ly] = in[gy * width + gx];
    }
    if (threadIdx.y < radius) {
        tile[threadIdx.y] = (gy >= radius) ? in[(gy - radius) * width + gx] : 0.0f;
        tile[ly + blockDim.y] = (gy + blockDim.y < height) ?
            in[(gy + blockDim.y) * width + gx] : 0.0f;
    }

    __syncthreads();

    if (gy < height) {
        float sum = 0.0f;
        for (int k = -radius; k <= radius; k++) {
            sum += c_kernelCol[k + radius] * tile[ly + k];
        }
        out[gy * width + gx] = sum;
    }
}

6. 3D模板（7点拉普拉斯）

#define TILE_X 16
#define TILE_Y 16
#define TILE_Z 4

__global__ void laplacian3D(
    float* out, const float* in,
    int nx, int ny, int nz
) {
    __shared__ float current[TILE_Y + 2][TILE_X + 2];
    __shared__ float above[TILE_Y][TILE_X];
    __shared__ float below[TILE_Y][TILE_X];

    int gx = blockIdx.x * TILE_X + threadIdx.x;
    int gy = blockIdx.y * TILE_Y + threadIdx.y;
    int gz = blockIdx.z * TILE_Z;

    int lx = threadIdx.x + 1;
    int ly = threadIdx.y + 1;

    // 处理TILE_Z平面
    for (int z = gz; z < min(gz + TILE_Z, nz - 1); z++) {
        if (z == 0) continue;

        // 加载当前平面带光环
        if (gx < nx && gy < ny) {
            current[ly][lx] = in[z * ny * nx + gy * nx + gx];
        }

        // 加载光环
        if (threadIdx.x == 0 && gx > 0) {
            current[ly][0] = in[z * ny * nx + gy * nx + (gx - 1)];
        }
        if (threadIdx.x == TILE_X - 1 && gx < nx - 1) {
            current[ly][TILE_X + 1] = in[z * ny * nx + gy * nx + (gx + 1)];
        }
        if (threadIdx.y == 0 && gy > 0) {
            current[0][lx] = in[z * ny * nx + (gy - 1) * nx + gx];
        }
        if (threadIdx.y == TILE_Y - 1 && gy < ny - 1) {
            current[TILE_Y + 1][lx] = in[z * ny * nx + (gy + 1) * nx + gx];
        }

        // 加载上下平面
        if (gx < nx && gy < ny) {
            above[threadIdx.y][threadIdx.x] = in[(z + 1) * ny * nx + gy * nx + gx];
            below[threadIdx.y][threadIdx.x] = in[(z - 1) * ny * nx + gy * nx + gx];
        }

        __syncthreads();

        // 计算7点模板
        if (gx >= 1 && gx < nx - 1 && gy >= 1 && gy < ny - 1) {
            out[z * ny * nx + gy * nx + gx] =
                -6.0f * current[ly][lx]
                + current[ly][lx - 1]   // x-1
                + current[ly][lx + 1]   // x+1
                + current[ly - 1][lx]   // y-1
                + current[ly + 1][lx]   // y+1
                + above[threadIdx.y][threadIdx.x]    // z+1
                + below[threadIdx.y][threadIdx.x];   // z-1
        }

        __syncthreads();
    }
}

7. 时间阻塞（多时间步）

// 在写回全局内存之前处理多个时间步
template <int TIMESTEPS>
__global__ void stencil_temporal_blocking(
    float* out, const float* in,
    int width, int height
) {
    // 更大的共享内存以适应时间扩展
    // 每个时间步将光环扩展1
    const int HALO = TIMESTEPS;
    extern __shared__ float smem[];

    float* current = smem;
    float* next = smem + (blockDim.y + 2 * HALO) * (blockDim.x + 2 * HALO);

    // 加载初始数据带扩展光环
    // ...

    __syncthreads();

    // 共享内存中的多个时间步
    for (int t = 0; t < TIMESTEPS; t++) {
        int shrinkHalo = TIMESTEPS - t - 1;
        int validXStart = shrinkHalo;
        int validXEnd = blockDim.x + 2 * HALO - shrinkHalo;
        int validYStart = shrinkHalo;
        int validYEnd = blockDim.y + 2 * HALO - shrinkHalo;

        int lx = threadIdx.x + HALO;
        int ly = threadIdx.y + HALO;

        // 只有在有效区域的线程计算
        if (lx >= validXStart + 1 && lx < validXEnd - 1 &&
            ly >= validYStart + 1 && ly < validYEnd - 1) {

            int PITCH = blockDim.x + 2 * HALO;
            next[ly * PITCH + lx] = -4.0f * current[ly * PITCH + lx]
                                  + current[ly * PITCH + lx - 1]
                                  + current[ly * PITCH + lx + 1]
                                  + current[(ly - 1) * PITCH + lx]
                                  + current[(ly + 1) * PITCH + lx];
        }

        __syncthreads();

        // 交换缓冲区
        float* temp = current;
        current = next;
        next = temp;

        __syncthreads();
    }

    // 将最终结果写入全局内存
    // ...
}

8. 边界条件模式

// 不同的边界条件策略
enum BoundaryCondition {
    BC_ZERO,       // 零填充
    BC_REPLICATE,  // 复制边缘值
    BC_REFLECT,    // 镜像反射
    BC_PERIODIC    // 环绕
};

__device__ inline int applyBoundary(int idx, int size, BoundaryCondition bc) {
    if (idx >= 0 && idx < size) return idx;

    switch (bc) {
        case BC_ZERO:
            return -1;  // 信号使用零
        case BC_REPLICATE:
            return (idx < 0) ? 0 : size - 1;
        case BC_REFLECT:
            if (idx < 0) return -idx - 1;
            if (idx >= size) return 2 * size - idx - 1;
            return idx;
        case BC_PERIODIC:
            return ((idx % size) + size) % size;
        default:
            return idx;
    }
}

__device__ inline float loadWithBoundary(
    const float* data, int x, int y,
    int width, int height, BoundaryCondition bc
) {
    int bx = applyBoundary(x, width, bc);
    int by = applyBoundary(y, height, bc);

    if (bx < 0 || by < 0) return 0.0f;

    return data[by * width + bx];
}

最佳实践

内存访问模式

模式	影响	建议
合并的全局读取	高	使线程访问与内存布局对齐
共享内存银行冲突	中	填充共享内存数组
光环加载效率	中	使用合作加载

平铺尺寸选择

GPU架构	推荐的平铺尺寸
Volta/Turing	32x32或16x16
Ampere	32x32
Hopper	32x32或64x32

性能提示

1. 使用常量内存 存储模板权重
2. 最大化共享内存中的数据重用
3. 考虑可分离的滤波器 对于2D卷积
4. 时间阻塞 对于迭代模板
5. 分析内存带宽 - 模板是内存受限的

流程集成

这项技能与以下流程集成：

• stencil-computation-optimization.js - 模板优化工作流
• gpu-image-video-processing.js - 图像过滤
• parallel-algorithm-design.js - 算法模式

输出格式

执行操作时，提供结构化输出：

{
  "operation": "generate-stencil",
  "status": "success",
  "stencil": {
    "type": "2D",
    "points": 5,
    "radius": 1,
    "boundary": "replicate"
  },
  "optimization": {
    "tile_size": [32, 32],
    "shared_memory_bytes": 4624,
    "halo_size": 1,
    "temporal_blocking": false
  },
  "performance": {
    "achieved_bandwidth_gbps": 850,
    "peak_bandwidth_gbps": 1555,
    "efficiency_percent": 54.7
  },
  "recommendations": [
    "Consider separable implementation for Gaussian filter",
    "Temporal blocking could reduce memory traffic by 2x"
  ],
  "artifacts": ["stencil_kernel.cu", "benchmark_results.json"]
}

约束

• 模板通常是内存带宽受限的
• 共享内存限制平铺尺寸
• 边界处理增加复杂性
• 3D模板具有更高的内存需求
• 分析以确保内存合并