The AI CUDA Engineer 👷

import torch
import torch.nn as nn
import torch.nn.functional as F


def module_fn(
    x: torch.Tensor,
    stride: int,
    padding: int,
    conv_weight: torch.Tensor,
    conv_bias: torch.Tensor,
) -> torch.Tensor:
    """
    Applies 3D convolution, max pooling, log sum exp, and ReLU activation.

    Args:
        x (torch.Tensor): Input tensor of shape (batch_size, in_channels, depth, height, width)
        stride (int): Stride of the convolution
        padding (int): Padding of the convolution
        conv_weight (torch.Tensor): Convolution weight tensor
        conv_bias (torch.Tensor): Convolution bias tensor

    Returns:
        torch.Tensor: Output tensor after applying convolution, max pooling, logsumexp and ReLU
    """
    x = F.conv3d(x, conv_weight, bias=conv_bias, stride=stride, padding=padding)
    x = F.max_pool3d(x, kernel_size=2, stride=2)
    x = torch.logsumexp(x, dim=1, keepdim=True)
    x = F.relu(x)
    return x


class Model(nn.Module):
    """
    Model that performs a 3D convolution, max pooling, log sum exp, and ReLU activation.
    """

    def __init__(self, in_channels, out_channels, kernel_size, stride, padding):
        super(Model, self).__init__()
        conv = nn.Conv3d(
            in_channels, out_channels, kernel_size, stride=stride, padding=padding
        )
        self.conv_weight = nn.Parameter(conv.weight)
        self.conv_bias = nn.Parameter(
            conv.bias
            + torch.randn(
                conv.bias.shape, device=conv.bias.device, dtype=conv.bias.dtype
            )
            * 0.02
        )

    def forward(self, x, stride, padding, fn=module_fn):
        return fn(x, stride, padding, self.conv_weight, self.conv_bias)


batch_size = 128
in_channels = 3
out_channels = 16
depth, height, width = 16, 32, 32
kernel_size = 3
stride = 1
padding = 1


def get_inputs():
    return [torch.randn(batch_size, in_channels, depth, height, width), stride, padding]


def get_init_inputs():
    return [in_channels, out_channels, kernel_size, stride, padding]

#include <torch/extension.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <cfloat>

template<int WARP_SIZE = 32>
__global__ void warp_uniform_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    const int N, const int C, const int D, const int H, const int W) {
    
    const int tid = threadIdx.x;
    const int wid = tid / WARP_SIZE;
    const int lane = tid % WARP_SIZE;
    const int bid = blockIdx.x;
    
    const int spatial_size = D * H * W;
    const int total_elements = N * spatial_size;
    const int warps_per_block = blockDim.x / WARP_SIZE;
    const int elements_per_iter = gridDim.x * warps_per_block * WARP_SIZE;
    
    // Process elements in warp-sized chunks
    for (int base_idx = bid * warps_per_block * WARP_SIZE + wid * WARP_SIZE + lane;
         base_idx < total_elements;
         base_idx += elements_per_iter) {
        
        // Calculate position (guaranteed uniform across warp)
        const int n = base_idx / spatial_size;
        const int pos = base_idx % spatial_size;
        const int d = pos / (H * W);
        const int hw = pos % (H * W);
        const int h = hw / W;
        const int w = hw % W;
        
        // Calculate base offset for this position
        const int base_offset = n * (C * spatial_size) + d * (H * W) + h * W + w;
        
        // Find maximum value (uniform across warp)
        float max_val = -FLT_MAX;
        
        #pragma unroll 4
        for (int c = 0; c < C; c++) {
            const float val = input[base_offset + c * spatial_size];
            max_val = fmaxf(max_val, val);
        }
        
        // Compute sum of exponentials (uniform across warp)
        float sum_exp = 0.0f;
        
        #pragma unroll 4
        for (int c = 0; c < C; c++) {
            const float val = input[base_offset + c * spatial_size];
            sum_exp += __expf(val - max_val);
        }
        
        // Compute final result (uniform across warp)
        const float result = fmaxf(0.0f, max_val + __logf(sum_exp));
        
        // Write result (coalesced write)
        if (base_idx < total_elements) {
            output[base_idx] = result;
        }
    }
}

torch::Tensor forward(
    torch::Tensor x,
    int64_t stride,
    int64_t padding,
    torch::Tensor conv_weight,
    torch::Tensor conv_bias) {

    x = x.contiguous();
    conv_weight = conv_weight.contiguous();
    conv_bias = conv_bias.contiguous();

    auto conv_result = torch::conv3d(x, conv_weight, conv_bias,
                                   {stride, stride, stride},
                                   {padding, padding, padding});
    
    auto pool_result = torch::max_pool3d(conv_result, {2, 2, 2}, {2, 2, 2});

    const int N = pool_result.size(0);
    const int C = pool_result.size(1);
    const int D = pool_result.size(2);
    const int H = pool_result.size(3);
    const int W = pool_result.size(4);

    auto output = torch::empty({N, 1, D, H, W}, pool_result.options());

    // Configure kernel launch parameters
    const int threads_per_block = 256; // Multiple of warp size (32)
    const int num_blocks = (N * D * H * W + threads_per_block - 1) / threads_per_block;
    
    warp_uniform_kernel<32><<<num_blocks, threads_per_block>>>(
        pool_result.data_ptr<float>(),
        output.data_ptr<float>(),
        N, C, D, H, W
    );

    return output;
}

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.def("forward", &forward, "Warp uniform kernel implementation");
}