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 <vector>
#include <cfloat>
#include <cmath>

__global__ void coalesced_memory_access_fused_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    const int N, const int C, const int D, const int H, const int W) {
    
    // Calculate total elements and stride for thread processing
    const int total_elements = N * D * H * W;
    const int num_threads = blockDim.x * gridDim.x;
    const int thread_id = blockIdx.x * blockDim.x + threadIdx.x;
    const int stride = D * H * W;
    
    // Each thread processes multiple elements using stride loop
    for (int idx = thread_id; idx < total_elements; idx += num_threads) {
        // Decode indices
        const int w = idx % W;
        int temp = idx / W;
        const int h = temp % H;
        temp /= H;
        const int d = temp % D;
        const int n = temp / D;

        // Calculate base index for current position
        const int base_idx = n * (C * stride) + d * (H * W) + h * W + w;

        // First pass: find maximum value across channels
        float max_val = -FLT_MAX;
        #pragma unroll
        for (int c = 0; c < C; ++c) {
            const int input_idx = base_idx + c * stride;
            max_val = fmaxf(max_val, input[input_idx]);
        }

        // Second pass: compute sum of exponentials
        float sum_exp = 0.0f;
        #pragma unroll
        for (int c = 0; c < C; ++c) {
            const int input_idx = base_idx + c * stride;
            sum_exp += expf(input[input_idx] - max_val);
        }

        // Compute final result with ReLU
        float result = max_val + logf(sum_exp);
        result = fmaxf(0.0f, result);
        
        // Write to output
        output[idx] = result;
    }
}

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

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

    // Perform max pooling using PyTorch
    auto pool_result = torch::max_pool3d(conv_result, {2, 2, 2}, {2, 2, 2});

    // Get dimensions for the fused logsumexp and ReLU operations
    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);

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

    // Launch kernel with stride-based processing
    const int block_size = 256;
    const int num_blocks = (N * D * H * W + block_size - 1) / block_size;
    
    coalesced_memory_access_fused_kernel<<<num_blocks, block_size>>>(
        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, "Coalesced memory access fused kernel");
}