```python
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.cpp_extension import load_inline


_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor relu_group_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>

constexpr int THREADS = 256;
constexpr int CHUNK = 4096;

__device__ __forceinline__ float relu_value(float x) {
    return x < 0.0f ? 0.0f : x;
}

__device__ __forceinline__ void block_sum(
    float& sum, float& square_sum
) {
    __shared__ float sums[8];
    __shared__ float square_sums[8];

    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        sum += __shfl_down_sync(0xffffffff, sum, offset);
        square_sum += __shfl_down_sync(0xffffffff, square_sum, offset);
    }

    if (lane == 0) {
        sums[warp] = sum;
        square_sums[warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = lane < 8 ? sums[lane] : 0.0f;
        square_sum = lane < 8 ? square_sums[lane] : 0.0f;

        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1) {
            sum += __shfl_down_sync(0xffffffff, sum, offset);
            square_sum += __shfl_down_sync(
                0xffffffff, square_sum, offset);
        }
    }
}

__global__ void partial_stats_kernel(
    const float* __restrict__ x,
    float2* __restrict__ partials,
    int64_t group_size,
    int chunks
) {
    const int group = blockIdx.x / chunks;
    const int chunk = blockIdx.x % chunks;
    const int64_t start = static_cast<int64_t>(chunk) * CHUNK;
    const int64_t base = static_cast<int64_t>(group) * group_size;

    float sum = 0.0f;
    float square_sum = 0.0f;

    #pragma unroll
    for (int k = 0; k < CHUNK / THREADS; ++k) {
        const int64_t i = start + threadIdx.x + k * THREADS;
        if (i < group_size) {
            const float value = relu_value(x[base + i]);
            sum += value;
            square_sum = fmaf(value, value, square_sum);
        }
    }

    block_sum(sum, square_sum);
    if (threadIdx.x == 0) {
        partials[blockIdx.x] = make_float2(sum, square_sum);
    }
}

__global__ void finish_stats_kernel(
    const float2* __restrict__ partials,
    float2* __restrict__ stats,
    int64_t group_size,
    int chunks,
    float eps
) {
    const int group = blockIdx.x;
    float sum = 0.0f;
    float square_sum = 0.0f;

    for (int i = threadIdx.x; i < chunks; i += THREADS) {
        const float2 value =
            partials[static_cast<int64_t>(group) * chunks + i];
        sum += value.x;
        square_sum += value.y;
    }

    block_sum(sum, square_sum);
    if (threadIdx.x == 0) {
        const float reciprocal = 1.0f / static_cast<float>(group_size);
        const float mean = sum * reciprocal;
        const float variance =
            fmaxf(fmaf(-mean, mean, square_sum * reciprocal), 0.0f);
        stats[group] = make_float2(mean, rsqrtf(variance + eps));
    }
}

__global__ void normalize_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float2* __restrict__ stats,
    float* __restrict__ output,
    int64_t count,
    int64_t spatial_size,
    int64_t group_size,
    int channels
) {
    const int64_t start =
        static_cast<int64_t>(blockIdx.x) * THREADS * 4 + threadIdx.x;

    #pragma unroll
    for (int k = 0; k < 4; ++k) {
        const int64_t i = start + k * THREADS;
        if (i < count) {
            const int group = i / group_size;
            const int channel = (i / spatial_size) % channels;
            const float2 parameters = stats[group];
            const float normalized =
                (relu_value(x[i]) - parameters.x) * parameters.y;
            output[i] = normalized * weight[channel] + bias[channel];
        }
    }
}

torch::Tensor relu_group_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps
) {
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");
    TORCH_CHECK(groups > 0 && x.size(1) % groups == 0,
                "Invalid group count");
    TORCH_CHECK(weight.device() == x.device() &&
                bias.device() == x.device() &&
                weight.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat &&
                weight.is_contiguous() && bias.is_contiguous() &&
                weight.numel() == x.size(1) &&
                bias.numel() == x.size(1),
                "Invalid affine parameters");

    const c10::cuda::CUDAGuard guard(x.device());
    auto output = torch::empty_like(x);
    if (x.numel() == 0) {
        return output;
    }

    const int channels = static_cast<int>(x.size(1));
    const int64_t spatial_size = x.size(2) * x.size(3) * x.size(4);
    const int64_t group_size = (channels / groups) * spatial_size;
    const int chunks = static_cast<int>(
        (group_size + CHUNK - 1) / CHUNK);
    const int64_t total_groups = x.size(0) * groups;

    auto partials = torch::empty({total_groups, chunks, 2}, x.options());
    auto stats = torch::empty({total_groups, 2}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    partial_stats_kernel<<<total_groups * chunks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float2*>(partials.data_ptr<float>()),
        group_size,
        chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_stats_kernel<<<total_groups, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partials.data_ptr<float>()),
        reinterpret_cast<float2*>(stats.data_ptr<float>()),
        group_size,
        chunks,
        static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const int64_t blocks =
        (x.numel() + THREADS * 4 - 1) / (THREADS * 4);
    normalize_kernel<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        reinterpret_cast<const float2*>(stats.data_ptr<float>()),
        output.data_ptr<float>(),
        x.numel(),
        spatial_size,
        group_size,
        channels);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_fused_ops = load_inline(
    name="fused_relu_group_norm_3d",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["relu_group_norm_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_channels, out_channels, kernel_size, groups, bias=False
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size, bias=bias
        )
        self.relu = nn.ReLU()
        self.group_norm = nn.GroupNorm(groups, out_channels)

    def forward(self, x):
        x = self.conv_transpose(x)
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
        ):
            return _fused_ops.relu_group_norm_cuda(
                x,
                self.group_norm.weight,
                self.group_norm.bias,
                self.group_norm.num_groups,
                self.group_norm.eps,
            )
        return self.group_norm(F.relu(x))
```