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]

/*
   This kernel fuses the logsumexp reduction over the channel dimension with ReLU activation.
   It uses a grid-stride loop to handle all output elements, and applies loop unrolling
   for the inner channel loop to improve instruction-level parallelism and reduce loop overhead.
*/

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

// Kernel definition
__global__ void fused_gridstride_logsumexp_relu_kernel(const float* __restrict__ input,
                                                         float* __restrict__ output,
                                                         const int N, const int C,
                                                         const int D, const int H, const int W) {
    // Total number of output elements (each corresponds to a reduction over channels)
    const int total = N * D * H * W;
    // spatial_stride is the stride between different channels for a fixed (n,d,h,w)
    const int spatial_stride = D * H * W;

    // Compute global thread id and total number of threads
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int num_threads = blockDim.x * gridDim.x;

    // Grid-stride loop: each thread processes multiple (n, d, h, w) positions
    for (int idx = tid; idx < total; idx += num_threads) {
        // Decode the linear index into (n, d, h, w): note that the tensor layout is [N, C, D, H, W]
        int tmp = idx;
        int w = tmp % W;
        tmp /= W;
        int h = tmp % H;
        tmp /= H;
        int d = tmp % D;
        int n = tmp / D;

        // Compute the base index for the reduction over channels
        int base = n * (C * spatial_stride) + d * (H * W) + h * W + w;

        // First pass: compute the maximum over channels for numerical stability
        float local_max = -FLT_MAX;
        #pragma unroll
        for (int c = 0; c < C; ++c) {
            float val = input[base + c * spatial_stride];
            local_max = fmaxf(local_max, val);
        }

        // Second pass: accumulate the sum of exponentials
        float sum_exp = 0.0f;
        #pragma unroll
        for (int c = 0; c < C; ++c) {
            float val = input[base + c * spatial_stride];
            sum_exp += expf(val - local_max);
        }

        // Compute the logsumexp result and apply ReLU
        float result = local_max + logf(sum_exp);
        result = fmaxf(0.0f, result);

        // Store the result into the output tensor
        output[idx] = result;
    }
}

// Host function: performs conv3d and max_pool3d using PyTorch, then fuses logsumexp and ReLU using the custom kernel

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

    // Ensure the tensors are contiguous
    x = x.contiguous();
    conv_weight = conv_weight.contiguous();
    conv_bias = conv_bias.contiguous();

    // Execute 3D convolution using PyTorch
    auto conv_result = torch::conv3d(x, conv_weight, conv_bias, {stride, stride, stride}, {padding, padding, padding});

    // Execute 3D max pooling using PyTorch with kernel size 2 and stride 2
    auto pool_result = torch::max_pool3d(conv_result, {2, 2, 2}, {2, 2, 2});

    // Dimensions of the pooled tensor are [N, C, D, H, W]
    int N = pool_result.size(0);
    int C = pool_result.size(1);
    int D = pool_result.size(2);
    int H = pool_result.size(3);
    int W = pool_result.size(4);

    // Create output tensor with shape [N, 1, D, H, W]
    auto output = torch::empty({N, 1, D, H, W}, pool_result.options());

    // Total number of output elements
    int total = N * D * H * W;
    const int block_size = 256;
    int num_blocks = (total + block_size - 1) / block_size;

    // Launch the fused kernel with grid-stride loop and unrolling for improved efficiency
    fused_gridstride_logsumexp_relu_kernel<<<num_blocks, block_size>>>(
        pool_result.data_ptr<float>(),
        output.data_ptr<float>(),
        N, C, D, H, W
    );
    cudaDeviceSynchronize();

    return output;
}

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.def("forward", &forward, "Fused grid-stride kernel for logsumexp and ReLU activation");
}