Agent skill

stencil-convolution

Expert skill for optimized stencil and convolution pattern implementations on GPU. Design tiled stencil algorithms with halos, implement 2D/3D convolution kernels, optimize boundary condition handling, apply temporal blocking techniques, generate separable filter implementations, and profile stencil memory bandwidth.

View SKILL.md on GitHub Repository
Stars 514
Forks 31

Install this agent skill to your Project

npx add-skill https://github.com/a5c-ai/babysitter/tree/main/library/specializations/gpu-programming/skills/stencil-convolution

Metadata

Additional technical details for this skill

author
babysitter-sdk
version
1.0.0
category
domain-algorithms
backlog id
SK-013

SKILL.md

stencil-convolution

You are stencil-convolution - a specialized skill for optimized stencil and convolution pattern implementations on GPU. This skill provides expert capabilities for scientific computing, image processing, and numerical simulations requiring neighborhood computations.

Overview

This skill enables AI-powered stencil and convolution operations including:

  • Designing tiled stencil algorithms with halos
  • Implementing 2D/3D convolution kernels
  • Optimizing boundary condition handling
  • Applying temporal blocking techniques
  • Generating separable filter implementations
  • Configuring shared memory tiling strategies
  • Profiling stencil memory bandwidth
  • Supporting multi-resolution stencils

Prerequisites

  • NVIDIA CUDA Toolkit 11.0+
  • GPU with compute capability 3.5+
  • Understanding of memory coalescing patterns
  • Nsight Compute for memory analysis
  • Optional: cuDNN for optimized convolutions

Capabilities

1. Basic 2D Stencil (5-Point Laplacian)

cuda
// Naive 5-point stencil (for comparison)
__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]      // left
                 + in[idx + 1]      // right
                 + in[idx - width]  // up
                 + in[idx + width]; // down
    }
}

2. Tiled Stencil with Shared Memory and Halo

cuda
#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 memory with halo
    __shared__ float tile[TILE_Y + 2 * HALO][TILE_X + 2 * HALO];

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

    // Local coordinates in shared memory (offset by halo)
    int lx = threadIdx.x + HALO;
    int ly = threadIdx.y + HALO;

    // Load center tile
    if (gx < width && gy < height) {
        tile[ly][lx] = in[gy * width + gx];
    }

    // Load halo regions
    // Left halo
    if (threadIdx.x < HALO && gx >= HALO) {
        tile[ly][lx - HALO] = in[gy * width + (gx - HALO)];
    }
    // Right halo
    if (threadIdx.x >= TILE_X - HALO && gx + HALO < width) {
        tile[ly][lx + HALO] = in[gy * width + (gx + HALO)];
    }
    // Top halo
    if (threadIdx.y < HALO && gy >= HALO) {
        tile[ly - HALO][lx] = in[(gy - HALO) * width + gx];
    }
    // Bottom halo
    if (threadIdx.y >= TILE_Y - HALO && gy + HALO < height) {
        tile[ly + HALO][lx] = in[(gy + HALO) * width + gx];
    }

    // Corner halos (if needed for larger stencils)
    // ...

    __syncthreads();

    // Compute stencil using shared memory
    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. Generic N-Point Stencil with Configurable Radius

cuda
template <int RADIUS>
__global__ void stencil2D_generic(
    float* out, const float* in,
    const float* weights,  // Stencil 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;

    // Load tile with halos
    // ... (similar to above but generic)

    __syncthreads();

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

        float result = 0.0f;
        int wIdx = 0;

        // Apply stencil weights
        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 Convolution with Arbitrary Kernel

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

// Kernel weights in constant memory for fast access
__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;

    // Load center
    if (gx < width && gy < height) {
        tile[ly * TILE_PITCH + lx] = in[gy * width + gx];
    } else {
        tile[ly * TILE_PITCH + lx] = 0.0f;  // Zero padding
    }

    // Load halos with boundary handling
    // Left
    if (threadIdx.x < kernelRadius) {
        int srcX = gx - kernelRadius;
        tile[ly * TILE_PITCH + (lx - kernelRadius)] =
            (srcX >= 0 && gy < height) ? in[gy * width + srcX] : 0.0f;
    }
    // Right
    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;
    }
    // Top and bottom (similar pattern)
    // ...

    __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. Separable Convolution (2-Pass for Performance)

cuda
// Separable convolution is faster: O(2*r) vs O(r^2)
// First pass: horizontal convolution
__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;

    // Load with halos
    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;
    }
}

// Second pass: vertical convolution
__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;

    // Load with halos
    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 Stencil (7-Point Laplacian)

cuda
#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;

    // Process TILE_Z planes
    for (int z = gz; z < min(gz + TILE_Z, nz - 1); z++) {
        if (z == 0) continue;

        // Load current plane with halos
        if (gx < nx && gy < ny) {
            current[ly][lx] = in[z * ny * nx + gy * nx + gx];
        }

        // Load halos
        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];
        }

        // Load above and below planes
        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();

        // Compute 7-point stencil
        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. Temporal Blocking (Multi-Timestep)

cuda
// Process multiple timesteps before writing back to global memory
template <int TIMESTEPS>
__global__ void stencil_temporal_blocking(
    float* out, const float* in,
    int width, int height
) {
    // Larger shared memory to accommodate temporal expansion
    // Each timestep expands the halo by 1
    const int HALO = TIMESTEPS;
    extern __shared__ float smem[];

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

    // Load initial data with expanded halo
    // ...

    __syncthreads();

    // Multiple timesteps in shared memory
    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;

        // Only threads in valid region compute
        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();

        // Swap buffers
        float* temp = current;
        current = next;
        next = temp;

        __syncthreads();
    }

    // Write final result to global memory
    // ...
}

8. Boundary Condition Patterns

cuda
// Different boundary condition strategies
enum BoundaryCondition {
    BC_ZERO,       // Zero padding
    BC_REPLICATE,  // Replicate edge values
    BC_REFLECT,    // Mirror reflection
    BC_PERIODIC    // Wrap around
};

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

    switch (bc) {
        case BC_ZERO:
            return -1;  // Signal to use zero
        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];
}

Best Practices

Memory Access Patterns

Pattern Impact Recommendation
Coalesced global reads High Align thread access to memory layout
Shared memory bank conflicts Medium Pad shared memory arrays
Halo loading efficiency Medium Use cooperative loading

Tile Size Selection

GPU Architecture Recommended Tile Size
Volta/Turing 32x32 or 16x16
Ampere 32x32
Hopper 32x32 or 64x32

Performance Tips

  1. Use constant memory for stencil weights
  2. Maximize data reuse in shared memory
  3. Consider separable filters for 2D convolutions
  4. Temporal blocking for iterative stencils
  5. Profile memory bandwidth - stencils are memory-bound

Process Integration

This skill integrates with the following processes:

  • stencil-computation-optimization.js - Stencil optimization workflows
  • gpu-image-video-processing.js - Image filtering
  • parallel-algorithm-design.js - Algorithm patterns

Output Format

When executing operations, provide structured output:

json
{
  "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"]
}

Constraints

  • Stencils are typically memory-bandwidth bound
  • Shared memory limits tile size
  • Boundary handling adds complexity
  • 3D stencils have higher memory requirements
  • Profile to ensure memory coalescing

Recommended Agent Skills

Expand your agent's capabilities with these related and highly-rated skills.

a5c-ai/babysitter

gsd-tools

Central utility skill for GSD operations. Provides config parsing, slug generation, timestamps, path operations, and orchestrates calls to other specialized skills. Acts as the unified entry point that the original gsd-tools.cjs provided via its lib/ modules (commands, config, core, init).

514 31
Explore
a5c-ai/babysitter

model-profile-resolution

Resolve model profile (quality/balanced/budget) at orchestration start and map agents to specific models. Enables cost/quality tradeoffs by selecting appropriate AI models for each agent role.

514 31
Explore
a5c-ai/babysitter

verification-suite

Plan structure validation, phase completeness checks, reference integrity verification, and artifact existence confirmation. Provides the structured verification layer ensuring GSD artifacts are well-formed and complete.

514 31
Explore
a5c-ai/babysitter

state-management

STATE.md reading, writing, and field-level updates. Provides cross-session state persistence via .planning/STATE.md with structured fields for current task, completed phases, blockers, decisions, and quick tasks.

514 31
Explore
a5c-ai/babysitter

git-integration

Git commit patterns, formats, and conventions for GSD methodology. Provides atomic commits per task, structured commit messages, planning file commits, branch management, and milestone tag operations.

514 31
Explore
a5c-ai/babysitter

frontmatter-parsing

YAML frontmatter parsing and manipulation for .planning/ documents. Provides read, write, update, query, and validation operations on frontmatter blocks in GSD markdown artifacts.

514 31
Explore

Didn't find tool you were looking for?

Be as detailed as possible for better results