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


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

torch::Tensor conv3d_cuda(
    torch::Tensor input,
    torch::Tensor weight,
    c10::optional<torch::Tensor> bias);
"""

_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>

namespace {

constexpr int TILE_M = 32;
constexpr int TILE_N = 128;
constexpr int TILE_K = 16;
constexpr int REDUCTION = 81;

__global__ void conv3d_3x3x3_kernel(
    const float* __restrict__ input,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int channels_out,
    int depth,
    int height,
    int width,
    int output_height,
    int output_width,
    int spatial_size)
{
    __shared__ float weights_tile[TILE_M][TILE_K];
    __shared__ float input_tile[TILE_K][TILE_N];

    const int tid = threadIdx.x;
    const int lane = tid & 31;
    const int channel_lane = tid >> 5;
    const int channel_start = blockIdx.y * TILE_M;
    const int spatial_start = blockIdx.x * TILE_N;

    const int64_t input_batch_stride =
        int64_t(3) * depth * height * width;
    const int64_t output_batch_stride =
        int64_t(channels_out) * spatial_size;

    input += int64_t(blockIdx.z) * input_batch_stride;
    output += int64_t(blockIdx.z) * output_batch_stride;

    // Every thread loads one spatial position for eight reduction indices.
    const int load_column = tid & (TILE_N - 1);
    const int position = spatial_start + load_column;
    const bool valid_position = position < spatial_size;

    const int output_plane = output_height * output_width;
    const int z = position / output_plane;
    const int remainder = position - z * output_plane;
    const int y = remainder / output_width;
    const int x = remainder - y * output_width;

    float accumulators[4][4];

    #pragma unroll
    for (int m = 0; m < 4; ++m) {
        #pragma unroll
        for (int n = 0; n < 4; ++n) {
            accumulators[m][n] = 0.0f;
        }
    }

    #pragma unroll 1
    for (int reduction_start = 0;
         reduction_start < REDUCTION;
         reduction_start += TILE_K) {

        #pragma unroll
        for (int i = tid; i < TILE_M * TILE_K; i += 256) {
            const int local_channel = i / TILE_K;
            const int local_k = i % TILE_K;
            const int channel = channel_start + local_channel;
            const int k = reduction_start + local_k;

            weights_tile[local_channel][local_k] =
                (channel < channels_out && k < REDUCTION)
                    ? weight[int64_t(channel) * REDUCTION + k]
                    : 0.0f;
        }

        #pragma unroll
        for (int local_k = tid / TILE_N;
             local_k < TILE_K;
             local_k += 2) {

            const int k = reduction_start + local_k;
            float value = 0.0f;

            if (valid_position && k < REDUCTION) {
                const int input_channel = k / 27;
                const int kernel_position = k - input_channel * 27;
                const int dz = kernel_position / 9;
                const int dy = (kernel_position % 9) / 3;
                const int dx = kernel_position % 3;

                const int64_t offset =
                    ((int64_t(input_channel) * depth + z + dz)
                        * height + y + dy)
                        * width + x + dx;

                value = input[offset];
            }

            input_tile[local_k][load_column] = value;
        }

        __syncthreads();

        #pragma unroll
        for (int k = 0; k < TILE_K; ++k) {
            float a[4];
            float b[4];

            #pragma unroll
            for (int m = 0; m < 4; ++m) {
                a[m] = weights_tile[channel_lane + m * 8][k];
            }

            #pragma unroll
            for (int n = 0; n < 4; ++n) {
                b[n] = input_tile[k][lane + n * 32];
            }

            #pragma unroll
            for (int m = 0; m < 4; ++m) {
                #pragma unroll
                for (int n = 0; n < 4; ++n) {
                    accumulators[m][n] =
                        fmaf(a[m], b[n], accumulators[m][n]);
                }
            }
        }

        __syncthreads();
    }

    #pragma unroll
    for (int m = 0; m < 4; ++m) {
        const int channel = channel_start + channel_lane + m * 8;

        if (channel < channels_out) {
            const float bias_value = bias ? bias[channel] : 0.0f;

            #pragma unroll
            for (int n = 0; n < 4; ++n) {
                const int position_out = spatial_start + lane + n * 32;
                if (position_out < spatial_size) {
                    output[int64_t(channel) * spatial_size + position_out] =
                        accumulators[m][n] + bias_value;
                }
            }
        }
    }
}

} // namespace

torch::Tensor conv3d_cuda(
    torch::Tensor input,
    torch::Tensor weight,
    c10::optional<torch::Tensor> bias)
{
    TORCH_CHECK(input.is_cuda() && weight.is_cuda(),
                "Input and weight must be CUDA tensors");
    TORCH_CHECK(input.device() == weight.device(),
                "Input and weight must use the same device");
    TORCH_CHECK(input.scalar_type() == at::kFloat &&
                weight.scalar_type() == at::kFloat,
                "Expected float32 tensors");
    TORCH_CHECK(input.is_contiguous() && weight.is_contiguous(),
                "Expected contiguous tensors");
    TORCH_CHECK(input.dim() == 5 && input.size(1) == 3,
                "Expected NCDHW input with three channels");
    TORCH_CHECK(weight.dim() == 5 && weight.size(1) == 3 &&
                weight.size(2) == 3 && weight.size(3) == 3 &&
                weight.size(4) == 3,
                "Expected a 3x3x3 convolution kernel");
    TORCH_CHECK(input.size(2) >= 3 &&
                input.size(3) >= 3 &&
                input.size(4) >= 3,
                "Input spatial dimensions must be at least three");

    const c10::cuda::CUDAGuard device_guard(input.device());

    const int batch = input.size(0);
    const int channels_out = weight.size(0);
    const int depth = input.size(2);
    const int height = input.size(3);
    const int width = input.size(4);
    const int output_depth = depth - 2;
    const int output_height = height - 2;
    const int output_width = width - 2;
    const int64_t spatial64 =
        int64_t(output_depth) * output_height * output_width;

    TORCH_CHECK(spatial64 <= 2147483647LL,
                "Output spatial size exceeds the supported range");
    TORCH_CHECK(batch <= 65535 &&
                (int64_t(channels_out) + TILE_M - 1) / TILE_M <= 65535,
                "Convolution exceeds CUDA grid limits");

    const float* bias_pointer = nullptr;
    if (bias.has_value()) {
        const auto& b = bias.value();
        TORCH_CHECK(b.device() == input.device() &&
                    b.scalar_type() == at::kFloat &&
                    b.is_contiguous() &&
                    b.dim() == 1 &&
                    b.numel() == channels_out,
                    "Invalid bias tensor");
        bias_pointer = b.data_ptr<float>();
    }

    auto output = torch::empty(
        {batch, channels_out, output_depth, output_height, output_width},
        input.options());

    if (output.numel() == 0) {
        return output;
    }

    const int spatial_size = static_cast<int>(spatial64);
    const dim3 grid(
        (spatial64 + TILE_N - 1) / TILE_N,
        (channels_out + TILE_M - 1) / TILE_M,
        batch);

    conv3d_3x3x3_kernel<<<
        grid, 256, 0, at::cuda::getCurrentCUDAStream()>>>(
            input.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias_pointer,
            output.data_ptr<float>(),
            channels_out,
            depth,
            height,
            width,
            output_height,
            output_width,
            spatial_size);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = None


def _get_extension():
    global _extension
    if _extension is None:
        _extension = load_inline(
            name="conv3d_implicit_gemm_3channel_v1",
            cpp_sources=_cpp_source,
            cuda_sources=_cuda_source,
            functions=["conv3d_cuda"],
            extra_cflags=["-O3"],
            extra_cuda_cflags=["-O3"],
            verbose=False,
        )
    return _extension


class ModelNew(nn.Module):
    def __init__(
        self,
        in_channels: int,
        out_channels: int,
        kernel_size: int,
        stride: int = 1,
        padding: int = 0,
        dilation: int = 1,
        groups: int = 1,
        bias: bool = False,
    ):
        super().__init__()
        self.conv3d = nn.Conv3d(
            in_channels,
            out_channels,
            (kernel_size, kernel_size, kernel_size),
            stride=stride,
            padding=padding,
            dilation=dilation,
            groups=groups,
            bias=bias,
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        conv = self.conv3d
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and conv.weight.dtype == torch.float32
            and x.dim() == 5
            and x.shape[1] == 3
            and conv.in_channels == 3
            and conv.kernel_size == (3, 3, 3)
            and conv.stride == (1, 1, 1)
            and conv.padding == (0, 0, 0)
            and conv.dilation == (1, 1, 1)
            and conv.groups == 1
            and not torch.is_autocast_enabled()
        ):
            return _get_extension().conv3d_cuda(
                x.contiguous(),
                conv.weight.contiguous(),
                conv.bias,
            )
        return conv(x)
```