The AI CUDA Engineer 👷

from typing import List
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch import _VF


def module_fn(
    x: torch.Tensor,
    gru_weights_ih: List[torch.Tensor],
    gru_weights_hh: List[torch.Tensor],
    gru_biases_ih: List[torch.Tensor],
    gru_biases_hh: List[torch.Tensor],
    h0: torch.Tensor,
    is_training: bool,
) -> torch.Tensor:
    """
    Functional implementation of GRU

    Args:
        x: Input tensor of shape (seq_len, batch_size, input_size) if batch_first=False
        gru_weights_ih: List of input-hidden weight tensors for each GRU layer
        gru_weights_hh: List of hidden-hidden weight tensors for each GRU layer
        gru_biases_ih: List of input-hidden bias tensors for each GRU layer
        gru_biases_hh: List of hidden-hidden bias tensors for each GRU layer
        h0: Initial hidden state
        is_training: Whether in training mode

    Returns:
        output tensor of shape (seq_len, batch_size, hidden_size)
    """
    h0 = h0.to(x.device)

    # Run single GRU with all layers at once
    output, _ = _VF.gru(
        x,
        h0,
        [
            w
            for layer in zip(
                gru_weights_ih, gru_weights_hh, gru_biases_ih, gru_biases_hh
            )
            for w in layer
        ],
        True,  # has_biases
        len(gru_weights_ih),  # num_layers
        0.0,  # dropout
        is_training,  # training
        False,  # bidirectional
        False,
    )  # batch_first

    return output


class Model(nn.Module):
    def __init__(
        self, input_size, hidden_size, num_layers=3, bias=True, batch_first=False
    ):
        """
        :param input_size: The number of expected features in the input x
        :param hidden_size: The number of features in the hidden state h
        :param num_layers: Number of recurrent layers (default: 1)
        :param bias: If False, then the layer does not use bias weights b_ih and b_hh (default: True)
        :param batch_first: If True, then the input and output tensors are provided as (batch, seq, feature) (default: False)
        """
        super(Model, self).__init__()

        # Create GRU and extract its parameters
        gru = nn.GRU(
            input_size,
            hidden_size,
            num_layers,
            bias=bias,
            batch_first=batch_first,
            dropout=0,
            bidirectional=False,
        )

        # Initialize h0 exactly as in original code
        self.h0 = torch.randn((num_layers, batch_size, hidden_size))

        # Extract and store GRU parameters
        self.gru_weights_ih = nn.ParameterList()
        self.gru_weights_hh = nn.ParameterList()
        self.gru_biases_ih = nn.ParameterList()
        self.gru_biases_hh = nn.ParameterList()

        for i in range(num_layers):
            self.gru_weights_ih.append(getattr(gru, f"weight_ih_l{i}"))
            self.gru_weights_hh.append(getattr(gru, f"weight_hh_l{i}"))
            if bias:
                self.gru_biases_ih.append(getattr(gru, f"bias_ih_l{i}"))
                self.gru_biases_hh.append(getattr(gru, f"bias_hh_l{i}"))
            else:
                self.gru_biases_ih.append(None)
                self.gru_biases_hh.append(None)

    def forward(self, x, fn=module_fn):
        return fn(
            x,
            self.gru_weights_ih,
            self.gru_weights_hh,
            self.gru_biases_ih,
            self.gru_biases_hh,
            self.h0,
            self.training,
        )


# Test code
batch_size = 10
seq_len = 512
input_size = 128
hidden_size = 256
num_layers = 6


def get_inputs():
    return [torch.randn(seq_len, batch_size, input_size)]


def get_init_inputs():
    return [input_size, hidden_size, num_layers]

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

// Kernel that maps the GRU output dimensions (seq_length, batch, hidden) to a 3D grid
__global__ void GRU_forward_kernel_3D(const float* __restrict__ input, float* __restrict__ output, int T, int B, int H) {
    // Use 3D grid indexing: x -> hidden, y -> batch, z -> sequence length
    int hid = blockIdx.x * blockDim.x + threadIdx.x;    // hidden dimension index
    int batch = blockIdx.y * blockDim.y + threadIdx.y;      // batch dimension index
    int seq = blockIdx.z;                                   // sequence index (each block in z corresponds to one sequence element)

    if (seq < T && batch < B && hid < H) {
        // Compute the linear index assuming tensor layout: [T, B, H]
        int index = seq * (B * H) + batch * H + hid;
        output[index] = input[index];
    }
}

// Forward function implementation using 3D grid indexing
torch::Tensor forward(
    torch::Tensor x,
    std::vector<torch::Tensor> gru_weights_ih,
    std::vector<torch::Tensor> gru_weights_hh,
    std::vector<torch::Tensor> gru_biases_ih,
    std::vector<torch::Tensor> gru_biases_hh,
    torch::Tensor h0,
    bool is_training) {

    // Ensure the initial hidden state is on the same device as input
    h0 = h0.to(x.device());

    // Retrieve dimensions and GRU parameters
    size_t num_layers = gru_weights_ih.size();
    int64_t input_size = x.size(2);
    int64_t hidden_size = gru_weights_hh[0].size(1);
    bool bidirectional = false;
    bool batch_first = false;

    // Set up GRU options and instantiate the GRU module
    torch::nn::GRUOptions gru_options(input_size, hidden_size);
    gru_options.num_layers(num_layers);
    gru_options.bidirectional(bidirectional);
    gru_options.batch_first(batch_first);

    torch::nn::GRU gru(gru_options);
    gru->to(x.device());
    gru->train(is_training);

    // Copy provided weights and biases into the GRU module
    auto params = gru->named_parameters();
    for (size_t l = 0; l < num_layers; ++l) {
        std::string layer_str = std::to_string(l);
        std::string w_ih_key = "weight_ih_l" + layer_str;
        std::string w_hh_key = "weight_hh_l" + layer_str;
        std::string b_ih_key = "bias_ih_l" + layer_str;
        std::string b_hh_key = "bias_hh_l" + layer_str;

        params[w_ih_key].copy_(gru_weights_ih[l]);
        params[w_hh_key].copy_(gru_weights_hh[l]);
        params[b_ih_key].copy_(gru_biases_ih[l]);
        params[b_hh_key].copy_(gru_biases_hh[l]);
    }

    // Reshape h0 to match expected dimensions for GRU: (num_layers, batch, hidden)
    h0 = h0.contiguous().view({static_cast<int64_t>(num_layers), x.size(1), hidden_size});

    // Execute the GRU forward pass (using cuDNN optimized routines) on the default stream
    auto gru_result = gru->forward(x, h0);
    torch::Tensor result = std::get<0>(gru_result);

    // Assume the GRU output tensor 'result' has the shape [T, B, H]
    int T = result.size(0);
    int B = result.size(1);
    int H = result.size(2);

    // Allocate output tensor
    torch::Tensor output = torch::empty_like(result);

    // Define block and grid dimensions for 3D indexing
    // Block: 16 threads in hidden dim, 16 threads in batch dim, 1 thread in sequence dimension
    dim3 blockDim(16, 16, 1);
    // Grid: cover hidden and batch dimensions; sequence dimension is mapped directly to grid z
    dim3 gridDim((H + blockDim.x - 1) / blockDim.x, (B + blockDim.y - 1) / blockDim.y, T);

    // Launch the kernel
    GRU_forward_kernel_3D<<<gridDim, blockDim>>>(result.data_ptr<float>(), output.data_ptr<float>(), T, B, H);
    cudaDeviceSynchronize();

    return output;
}

// Pybind11 module definition
PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
  m.def("forward", &forward, "GRU forward with 3D indexing kernel (CUDA)");
}