The AI CUDA Engineer 👷

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


def module_fn(
    x: torch.Tensor,
    conv_transpose_weight: torch.Tensor,
    conv_transpose_bias: torch.Tensor,
    sum_weight: torch.Tensor,
    norm_weight: torch.Tensor,
    norm_bias: torch.Tensor,
    stride: tuple,
    padding: tuple,
    output_padding: tuple,
    pool_kernel_size: tuple,
    norm_shape: tuple,
) -> torch.Tensor:
    """
    Functional implementation of a sequence of operations:
    1. 3D transposed convolution
    2. Addition with a learnable weight
    3. Layer normalization
    4. 3D average pooling
    5. GELU activation

    Args:
        x (torch.Tensor): Input tensor of shape (batch_size, in_channels, depth, height, width)
        conv_transpose_weight (torch.Tensor): Weight tensor for transposed convolution
        conv_transpose_bias (torch.Tensor): Bias tensor for transposed convolution
        sum_weight (torch.Tensor): Learnable weight for addition
        norm_weight (torch.Tensor): Weight tensor for layer normalization
        norm_bias (torch.Tensor): Bias tensor for layer normalization
        stride (tuple): Stride for transposed convolution, as (depth_stride, height_stride, width_stride)
        padding (tuple): Padding for transposed convolution, as (depth_pad, height_pad, width_pad)
        output_padding (tuple): Output padding for transposed convolution, as (depth_pad, height_pad, width_pad)
        pool_kernel_size (tuple): Kernel size for average pooling, as (depth_kernel, height_kernel, width_kernel)
        norm_shape (tuple): Shape for layer normalization

    Returns:
        torch.Tensor: Output tensor after applying all operations
    """
    x = F.conv_transpose3d(
        x,
        conv_transpose_weight,
        bias=conv_transpose_bias,
        stride=stride,
        padding=padding,
        output_padding=output_padding,
    )
    x = x + sum_weight
    x = F.layer_norm(x, norm_shape, norm_weight, norm_bias)
    x = F.avg_pool3d(x, kernel_size=pool_kernel_size)
    x = F.gelu(x)
    return x


class Model(nn.Module):
    """
    Model that performs a 3D transposed convolution, followed by a sum, layer normalization, average pooling, and GELU activation.
    """

    def __init__(
        self,
        in_channels,
        out_channels,
        kernel_size,
        stride,
        padding,
        output_padding,
        sum_weight,
        norm_shape,
        pool_kernel_size,
    ):
        super(Model, self).__init__()
        conv = nn.ConvTranspose3d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
            output_padding=output_padding,
        )
        self.conv_transpose_weight = nn.Parameter(conv.weight)
        self.conv_transpose_bias = nn.Parameter(conv.bias)
        self.sum_weight = nn.Parameter(torch.tensor(sum_weight))
        norm = nn.LayerNorm(norm_shape)
        self.norm_weight = nn.Parameter(norm.weight + torch.randn(norm_shape) * 0.02)
        self.norm_bias = nn.Parameter(norm.bias + torch.randn(norm_shape) * 0.02)

    def forward(
        self,
        x,
        stride,
        padding,
        output_padding,
        pool_kernel_size,
        norm_shape,
        fn=module_fn,
    ):
        return fn(
            x,
            self.conv_transpose_weight,
            self.conv_transpose_bias,
            self.sum_weight,
            self.norm_weight,
            self.norm_bias,
            stride,
            padding,
            output_padding,
            pool_kernel_size,
            norm_shape,
        )


batch_size = 128
in_channels = 32
out_channels = 64
depth, height, width = 16, 32, 32
kernel_size = (3, 3, 3)
stride = (2, 2, 2)
padding = (1, 1, 1)
output_padding = (1, 1, 1)
sum_weight = 1.0
norm_shape = (out_channels,)
pool_kernel_size = (2, 2, 2)


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


def get_init_inputs():
    return [
        in_channels,
        out_channels,
        kernel_size,
        stride,
        padding,
        output_padding,
        sum_weight,
        norm_shape,
        pool_kernel_size,
    ]

#include <torch/extension.h>
#include <ATen/ATen.h>
#include <cuda_runtime.h>
#include <vector>

#define BLOCK_SIZE 256
#define WARP_SIZE 32

// This kernel implements layer normalization with manual loop unrolling
// to reduce loop overhead and improve performance. The critical loops that
// accumulate the sum and sum of squares, as well as the normalization loop,
// are unrolled with '#pragma unroll'.

template <typename T>
__global__ void unrolled_layernorm_kernel(
    const T* __restrict__ input,
    const T* __restrict__ gamma,
    const T* __restrict__ beta,
    T* __restrict__ output,
    int n1,
    int n2
) {
    // Shared memory for accumulating sum and square sum
    __shared__ float s_mean;
    __shared__ float s_variance;
    if (threadIdx.x == 0) {
        s_mean = 0.0f;
        s_variance = 0.0f;
    }
    __syncthreads();

    const int tid = threadIdx.x;
    const int bid = blockIdx.x;  // each block processes one normalization group
    const int offset = bid * n2;

    float local_sum = 0.0f;
    float local_sq_sum = 0.0f;

    // Manually unrolled loop for summing: assume n2 is small enough to benefit
    #pragma unroll
    for (int i = tid; i < n2; i += BLOCK_SIZE) {
        float val = __ldg(&input[offset + i]);
        local_sum += val;
        local_sq_sum += val * val;
    }

    // Warp-level reduction using shuffle instructions with unrolling
    float warp_sum = local_sum;
    float warp_sq_sum = local_sq_sum;
    #pragma unroll
    for (int offset_w = WARP_SIZE/2; offset_w > 0; offset_w /= 2) {
        warp_sum += __shfl_down_sync(0xffffffff, warp_sum, offset_w);
        warp_sq_sum += __shfl_down_sync(0xffffffff, warp_sq_sum, offset_w);
    }

    // Only the first thread of each warp performs atomic addition
    if ((tid & (WARP_SIZE - 1)) == 0) {
        atomicAdd(&s_mean, warp_sum);
        atomicAdd(&s_variance, warp_sq_sum);
    }
    __syncthreads();

    // Compute mean and variance; rsqrtf provides numerical stability
    float mean = s_mean / n2;
    float variance = s_variance / n2 - mean * mean;
    float inv_std = rsqrtf(variance + 1e-5f);
    
    // Unrolled loop for normalization and applying scale and bias
    #pragma unroll
    for (int i = tid; i < n2; i += BLOCK_SIZE) {
        int idx = offset + i;
        float normalized = (input[idx] - mean) * inv_std;
        output[idx] = gamma[i] * normalized + beta[i];
    }
}

// The forward function performs conv_transpose3d, adds sum_weight, applies our unrolled fused layer norm,
// then follows with avg_pool3d and GELU activation.

torch::Tensor forward(
    torch::Tensor x,
    torch::Tensor conv_transpose_weight,
    torch::Tensor conv_transpose_bias,
    torch::Tensor sum_weight,
    torch::Tensor norm_weight,
    torch::Tensor norm_bias,
    std::vector<int64_t> stride,
    std::vector<int64_t> padding,
    std::vector<int64_t> output_padding,
    std::vector<int64_t> pool_kernel_size,
    std::vector<int64_t> norm_shape
) {
    // Ensure all inputs are contiguous
    x = x.contiguous();
    conv_transpose_weight = conv_transpose_weight.contiguous();
    sum_weight = sum_weight.contiguous();
    norm_weight = norm_weight.contiguous();
    norm_bias = norm_bias.contiguous();

    at::IntArrayRef strideRef(stride);
    at::IntArrayRef paddingRef(padding);
    at::IntArrayRef outputPaddingRef(output_padding);
    at::IntArrayRef poolKernelRef(pool_kernel_size);

    // 1. 3D transposed convolution
    auto out = at::conv_transpose3d(
        x,
        conv_transpose_weight,
        conv_transpose_bias,
        strideRef,
        paddingRef,
        outputPaddingRef,
        /*groups=*/1,
        /*dilation=*/1
    );

    // 2. Elementwise addition with sum_weight
    out.add_(sum_weight);

    // 3. Custom fused layer normalization with unrolled loops
    auto out_size = out.sizes();
    int64_t n1 = 1;
    for (int i = 0; i < out.dim() - norm_shape.size(); ++i) {
        n1 *= out_size[i];
    }
    int64_t n2 = 1;
    for (size_t i = 0; i < norm_shape.size(); ++i) {
        n2 *= norm_shape[i];
    }

    auto output = torch::empty_like(out);
    dim3 grid(n1);
    dim3 block(BLOCK_SIZE);

    unrolled_layernorm_kernel<float><<<grid, block>>>(
        out.data_ptr<float>(),
        norm_weight.data_ptr<float>(),
        norm_bias.data_ptr<float>(),
        output.data_ptr<float>(),
        n1,
        n2
    );
    cudaDeviceSynchronize();

    // 4. 3D average pooling
    output = at::avg_pool3d(
        output,
        poolKernelRef,          // kernel_size
        poolKernelRef,          // stride (same as kernel_size)
        {0,0,0},
        false,
        true
    );

    // 5. GELU activation
    output = at::gelu(output);

    return output;
}

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.def("forward", &forward, "Unrolled fused layer norm forward (CUDA)");
}