The AI CUDA Engineer 👷

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


def module_fn(
    x: torch.Tensor,
    conv_weight: torch.Tensor,
    conv_bias: torch.Tensor,
) -> torch.Tensor:
    """
    Applies 3D convolution, HardSwish, ReLU, Softmax and mean reduction.

    Args:
        x (torch.Tensor): Input tensor of shape (batch_size, in_channels, depth, height, width)
        conv_weight (torch.Tensor): 3D convolution weight tensor of shape
            (out_channels, in_channels, kernel_size, kernel_size, kernel_size)
        conv_bias (torch.Tensor): Bias tensor for 3D convolution of shape (out_channels)

    Returns:
        torch.Tensor: Output tensor after applying convolution, activations and reduction,
            with shape (batch_size, out_channels)
    """
    x = F.conv3d(x, conv_weight, bias=conv_bias)
    x = F.hardswish(x)
    x = F.relu(x)
    x = F.softmax(x, dim=1)
    x = torch.mean(x, dim=[2, 3, 4])
    return x


class Model(nn.Module):
    """
    Simple model that performs a 3D convolution, applies HardSwish, ReLU, Softmax, and then calculates the mean.
    """

    def __init__(self, in_channels, out_channels, kernel_size):
        super(Model, self).__init__()

        conv = nn.Conv3d(in_channels, out_channels, kernel_size)
        self.conv_weight = nn.Parameter(conv.weight)
        self.conv_bias = nn.Parameter(conv.bias + torch.ones_like(conv.bias) * 0.02)

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


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


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


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

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

// Combined HardSwish and ReLU with aligned memory access
__global__ void fused_hswish_relu_kernel(const float* __restrict__ input,
                                        float* __restrict__ output,
                                        int64_t size) {
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index < size) {
        // Use __ldg for read-only access with 128-bit alignment assumption
        float x = __ldg(input + index);
        float relu6 = fminf(fmaxf(x + 3.0f, 0.0f), 6.0f);
        output[index] = fmaxf(x * relu6 * 0.16666667f, 0.0f);
    }
}

// Softmax kernel with read-only cache optimization
__global__ void optimized_softmax_kernel(const float* __restrict__ input,
                                      float* __restrict__ output,
                                      int batch_size,
                                      int channels,
                                      int spatial_size) {
    extern __shared__ float shared_max[];
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int tid = threadIdx.x;
    int total_elements = batch_size * spatial_size;

    if (index < total_elements) {
        int batch_idx = index / spatial_size;
        int spatial_idx = index % spatial_size;

        // Load through read-only cache with alignment
        float max_val = -FLT_MAX;
        for (int c = 0; c < channels; ++c) {
            int idx = batch_idx * channels * spatial_size + c * spatial_size + spatial_idx;
            max_val = fmaxf(max_val, __ldg(input + idx));
        }
        shared_max[tid] = max_val;

        // Block-wide max reduction
        __syncthreads();
        for (int s = blockDim.x/2; s > 0; s >>= 1) {
            if (tid < s) {
                shared_max[tid] = fmaxf(shared_max[tid], shared_max[tid + s]);
            }
            __syncthreads();
        }
        float block_max = shared_max[0];

        // Compute exp with stable values
        float sum_exp = 0.0f;
        for (int c = 0; c < channels; ++c) {
            int idx = batch_idx * channels * spatial_size + c * spatial_size + spatial_idx;
            float val = expf(__ldg(input + idx) - block_max);
            sum_exp += val;
            output[idx] = val;
        }

        // Normalize using fast reciprocal
        float inv_sum = 1.0f / sum_exp;
        for (int c = 0; c < channels; ++c) {
            int idx = batch_idx * channels * spatial_size + c * spatial_size + spatial_idx;
            output[idx] *= inv_sum;
        }
    }
}

torch::Tensor module_forward(
    torch::Tensor x,
    torch::Tensor conv_weight,
    torch::Tensor conv_bias)
{
    // Ensure 128-bit aligned memory access
    x = x.contiguous().cuda();
    conv_weight = conv_weight.contiguous().cuda();
    conv_bias = conv_bias.contiguous().cuda();

    x = torch::conv3d(x, conv_weight, conv_bias);

    int64_t batch_size = x.size(0);
    int64_t channels = x.size(1);
    int64_t depth = x.size(2);
    int64_t height = x.size(3);
    int64_t width = x.size(4);
    int64_t total_size = batch_size * channels * depth * height * width;

    torch::Tensor x_processed = torch::empty_like(x);

    // Launch fused HS-RReLU with optimal block size
    const int block_size = 256;
    int blocks = (total_size + block_size - 1) / block_size;
    fused_hswish_relu_kernel<<<blocks, block_size>>>(
        x.data_ptr<float>(),
        x_processed.data_ptr<float>(),
        total_size
    );

    // Process softmax
    x_processed = x_processed.view({batch_size, channels, depth * height * width});
    torch::Tensor x_softmax = torch::empty_like(x_processed);
    
    int spatial_size = depth * height * width;
    blocks = (batch_size * spatial_size + block_size - 1) / block_size;
    optimized_softmax_kernel<<<blocks, block_size, block_size * sizeof(float)>>>(
        x_processed.data_ptr<float>(),
        x_softmax.data_ptr<float>(),
        batch_size,
        channels,
        spatial_size
    );

    return x_softmax.view({batch_size, channels, depth, height, width})
                    .mean({2, 3, 4});
}

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.def("forward", &module_forward, "CUDA optimized memory forward");
}