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>

// Shared memory kernel for better performance
__global__ void optimized_fused_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    const int N, const int C, const int D, const int H, const int W) {
    
    extern __shared__ float shared_data[];
    
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int num_blocks = gridDim.x;
    const int stride = D * H * W;
    const int total_elements = N * D * H * W;
    
    // Process multiple elements per thread using grid-stride loop
    for (int idx = bid * blockDim.x + tid; idx < total_elements; idx += blockDim.x * num_blocks) {
        // 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;
        
        // Use shared memory for intermediate calculations
        float max_val = -FLT_MAX;
        float local_sum = 0.0f;
        
        // First pass: find maximum (coalesced memory access)
        #pragma unroll 4
        for (int c = 0; c < C; ++c) {
            const int input_idx = n * (C * stride) + c * stride + d * (H * W) + h * W + w;
            max_val = fmaxf(max_val, input[input_idx]);
        }
        
        // Store max_val in shared memory
        shared_data[tid] = max_val;
        __syncthreads();
        
        // Second pass: compute sum of exponentials
        #pragma unroll 4
        for (int c = 0; c < C; ++c) {
            const int input_idx = n * (C * stride) + c * stride + d * (H * W) + h * W + w;
            local_sum += __expf(input[input_idx] - shared_data[tid]);
        }
        
        // Compute final result with ReLU using intrinsics for better performance
        float result = shared_data[tid] + __logf(local_sum);
        result = fmaxf(0.0f, result);
        
        // Write to output (coalesced write)
        output[idx] = result;
    }
}

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

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

    // 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});

    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());

    // Optimize kernel launch configuration
    const int block_size = 256;
    const int num_blocks = std::min(65535, (N * D * H * W + block_size - 1) / block_size);
    
    optimized_fused_kernel<<<num_blocks, block_size, block_size * sizeof(float)>>>(
        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, "Optimized fused 3D operations");
}