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
#define WARPS_PER_BLOCK (BLOCK_SIZE/WARP_SIZE)

template <typename T>
__global__ void selective_atomic_layernorm_kernel(
    const T* __restrict__ input,
    const T* __restrict__ gamma,
    const T* __restrict__ beta,
    T* __restrict__ output,
    int n1,
    int n2
) {
    __shared__ float s_partial_sums[WARPS_PER_BLOCK];
    __shared__ float s_partial_squares[WARPS_PER_BLOCK];
    __shared__ float s_mean, s_variance;

    const int tid = threadIdx.x;
    const int wid = tid / WARP_SIZE;
    const int lane = tid % WARP_SIZE;
    const int bid = blockIdx.x;
    const int offset = bid * n2;

    // Initialize shared memory
    if (tid < WARPS_PER_BLOCK) {
        s_partial_sums[tid] = 0.0f;
        s_partial_squares[tid] = 0.0f;
    }
    if (tid == 0) {
        s_mean = 0.0f;
        s_variance = 0.0f;
    }
    __syncthreads();

    // Phase 1: Warp-level reduction without atomics
    float local_sum = 0.0f;
    float local_sq_sum = 0.0f;

    // Use vectorized loads where possible
    if (n2 >= 4) {
        const float4* input4 = reinterpret_cast<const float4*>(input + offset);
        for (int i = tid; i < n2/4; i += BLOCK_SIZE) {
            float4 vals = input4[i];
            local_sum += vals.x + vals.y + vals.z + vals.w;
            local_sq_sum += vals.x * vals.x + vals.y * vals.y + 
                           vals.z * vals.z + vals.w * vals.w;
        }
        
        // Handle remaining elements
        for (int i = (n2/4)*4 + tid; i < n2; i += BLOCK_SIZE) {
            float val = input[offset + i];
            local_sum += val;
            local_sq_sum += val * val;
        }
    } else {
        for (int i = tid; i < n2; i += BLOCK_SIZE) {
            float val = input[offset + i];
            local_sum += val;
            local_sq_sum += val * val;
        }
    }

    // Warp-level reduction using shuffle
    #pragma unroll
    for (int offset = WARP_SIZE/2; offset > 0; offset >>= 1) {
        local_sum += __shfl_down_sync(0xffffffff, local_sum, offset);
        local_sq_sum += __shfl_down_sync(0xffffffff, local_sq_sum, offset);
    }

    // First thread in each warp writes to shared memory
    if (lane == 0) {
        s_partial_sums[wid] = local_sum;
        s_partial_squares[wid] = local_sq_sum;
    }
    __syncthreads();

    // Phase 2: Final reduction using minimal atomics
    if (tid == 0) {
        float final_sum = 0.0f;
        float final_sq_sum = 0.0f;
        
        #pragma unroll
        for (int i = 0; i < WARPS_PER_BLOCK; ++i) {
            final_sum += s_partial_sums[i];
            final_sq_sum += s_partial_squares[i];
        }
        
        s_mean = final_sum / n2;
        s_variance = final_sq_sum / n2 - s_mean * s_mean;
    }
    __syncthreads();

    float mean = s_mean;
    float inv_std = rsqrtf(s_variance + 1e-5f);

    // Vectorized output computation where possible
    if (n2 >= 4) {
        float4* output4 = reinterpret_cast<float4*>(output + offset);
        const float4* input4 = reinterpret_cast<const float4*>(input + offset);
        const float4* gamma4 = reinterpret_cast<const float4*>(gamma);
        const float4* beta4 = reinterpret_cast<const float4*>(beta);
        
        for (int i = tid; i < n2/4; i += BLOCK_SIZE) {
            float4 in = input4[i];
            float4 g = gamma4[i];
            float4 b = beta4[i];
            
            float4 out;
            out.x = g.x * ((in.x - mean) * inv_std) + b.x;
            out.y = g.y * ((in.y - mean) * inv_std) + b.y;
            out.z = g.z * ((in.z - mean) * inv_std) + b.z;
            out.w = g.w * ((in.w - mean) * inv_std) + b.w;
            
            output4[i] = out;
        }
        
        // Handle remaining elements
        for (int i = (n2/4)*4 + tid; i < n2; i += BLOCK_SIZE) {
            float normalized = (input[offset + i] - mean) * inv_std;
            output[offset + i] = gamma[i] * normalized + beta[i];
        }
    } else {
        for (int i = tid; i < n2; i += BLOCK_SIZE) {
            float normalized = (input[offset + i] - mean) * inv_std;
            output[offset + i] = gamma[i] * normalized + beta[i];
        }
    }
}

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

    auto out = at::conv_transpose3d(
        x,
        conv_transpose_weight,
        conv_transpose_bias,
        strideRef,
        paddingRef,
        outputPaddingRef,
        1,
        1
    );

    out.add_(sum_weight);

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

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

    output = at::avg_pool3d(
        output,
        poolKernelRef,
        poolKernelRef,
        {0,0,0},
        false,
        true
    );

    output = at::gelu(output);

    return output;
}

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