```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 dilated_conv_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t sh, int64_t sw,
    int64_t ph, int64_t pw,
    int64_t dh, int64_t dw);
"""

_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>
#include <cublas_v2.h>
#include <algorithm>
#include <climits>

template<bool SPECIALIZED>
__global__ void unfold_kernel(
    const float* __restrict__ x,
    float* __restrict__ columns,
    int channels, int height, int width,
    int kh_arg, int kw_arg,
    int output_width,
    int sh_arg, int sw_arg,
    int ph, int pw,
    int dh_arg, int dw_arg,
    int first_position, int positions)
{
    const int kh = SPECIALIZED ? 5 : kh_arg;
    const int kw = SPECIALIZED ? 9 : kw_arg;
    const int sh = SPECIALIZED ? 1 : sh_arg;
    const int sw = SPECIALIZED ? 1 : sw_arg;
    const int dh = SPECIALIZED ? 2 : dh_arg;
    const int dw = SPECIALIZED ? 3 : dw_arg;

    const int reduction_index = blockIdx.y;
    const int channel = reduction_index / (kh * kw);
    const int kernel_position = reduction_index % (kh * kw);
    const int kernel_y = kernel_position / kw;
    const int kernel_x = kernel_position % kw;

    for (int p = blockIdx.x * blockDim.x + threadIdx.x;
         p < positions;
         p += blockDim.x * gridDim.x)
    {
        const int position = first_position + p;
        const int iy = (position / output_width) * sh - ph + kernel_y * dh;
        const int ix = (position % output_width) * sw - pw + kernel_x * dw;

        float value = 0.0f;
        if (static_cast<unsigned>(iy) < static_cast<unsigned>(height) &&
            static_cast<unsigned>(ix) < static_cast<unsigned>(width)) {
            value = x[
                (static_cast<int64_t>(channel) * height + iy) * width + ix
            ];
        }

        columns[static_cast<int64_t>(reduction_index) * positions + p] = value;
    }
}

__global__ void add_bias_kernel(
    float* __restrict__ output,
    const float* __restrict__ bias,
    int64_t count, int spatial, int channels)
{
    for (int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
         i < count;
         i += static_cast<int64_t>(blockDim.x) * gridDim.x)
    {
        output[i] += bias[(i / spatial) % channels];
    }
}

torch::Tensor dilated_conv_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t sh, int64_t sw,
    int64_t ph, int64_t pw,
    int64_t dh, int64_t dw)
{
    const c10::cuda::CUDAGuard device_guard(x.device());

    TORCH_CHECK(x.is_cuda() && weight.is_cuda(), "Expected CUDA tensors");
    TORCH_CHECK(x.scalar_type() == at::kFloat &&
                weight.scalar_type() == at::kFloat, "Expected float32 tensors");
    TORCH_CHECK(x.dim() == 4 && weight.dim() == 4, "Expected four dimensions");
    TORCH_CHECK(x.is_contiguous() && weight.is_contiguous(),
                "Expected contiguous tensors");
    TORCH_CHECK(x.device() == weight.device(), "Device mismatch");

    const int64_t batch = x.size(0);
    const int64_t channels = x.size(1);
    const int64_t height = x.size(2);
    const int64_t width = x.size(3);
    const int64_t out_channels = weight.size(0);
    const int64_t kh = weight.size(2);
    const int64_t kw = weight.size(3);

    TORCH_CHECK(weight.size(1) == channels, "Channel mismatch");
    TORCH_CHECK(sh > 0 && sw > 0 && dh > 0 && dw > 0, "Invalid convolution");
    TORCH_CHECK(ph >= 0 && pw >= 0, "Invalid padding");

    const int64_t effective_h = height + 2 * ph - dh * (kh - 1) - 1;
    const int64_t effective_w = width + 2 * pw - dw * (kw - 1) - 1;
    TORCH_CHECK(effective_h >= 0 && effective_w >= 0, "Kernel exceeds input");

    const int64_t oh = effective_h / sh + 1;
    const int64_t ow = effective_w / sw + 1;
    const int64_t spatial = oh * ow;
    const int64_t reduction = channels * kh * kw;

    TORCH_CHECK(spatial <= INT_MAX && out_channels <= INT_MAX &&
                reduction <= 65535 && height <= INT_MAX && width <= INT_MAX,
                "Convolution dimensions exceed supported limits");

    auto output = torch::empty({batch, out_channels, oh, ow}, x.options());
    if (output.numel() == 0) {
        return output;
    }

    constexpr int64_t workspace_elements = 48 * 1024 * 1024;
    int64_t tile_size = std::min<int64_t>(
        32768, std::max<int64_t>(256, workspace_elements / reduction));
    tile_size = std::max<int64_t>(256, (tile_size / 256) * 256);
    tile_size = std::min(tile_size, spatial);

    auto columns = torch::empty({reduction, tile_size}, x.options());

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    cublasHandle_t handle = at::cuda::getCurrentCUDABlasHandle();

    const float alpha = 1.0f;
    const float beta = 0.0f;

    const bool specialized =
        kh == 5 && kw == 9 &&
        sh == 1 && sw == 1 &&
        dh == 2 && dw == 3;

    for (int64_t n = 0; n < batch; ++n) {
        const float* input_ptr =
            x.data_ptr<float>() + n * channels * height * width;

        for (int64_t start = 0; start < spatial; start += tile_size) {
            const int positions =
                static_cast<int>(std::min(tile_size, spatial - start));

            const dim3 grid(
                std::min(32, (positions + 255) / 256),
                static_cast<unsigned>(reduction));

            if (specialized) {
                unfold_kernel<true><<<grid, 256, 0, stream>>>(
                    input_ptr, columns.data_ptr<float>(),
                    channels, height, width, kh, kw, ow,
                    sh, sw, ph, pw, dh, dw, start, positions);
            } else {
                unfold_kernel<false><<<grid, 256, 0, stream>>>(
                    input_ptr, columns.data_ptr<float>(),
                    channels, height, width, kh, kw, ow,
                    sh, sw, ph, pw, dh, dw, start, positions);
            }
            C10_CUDA_KERNEL_LAUNCH_CHECK();

            float* output_ptr =
                output.data_ptr<float>() + n * out_channels * spatial + start;

            // Interpret the row-major unfolded input and weights as
            // column-major matrices. Write directly into the NCHW output.
            const cublasStatus_t status = cublasGemmEx(
                handle,
                CUBLAS_OP_N, CUBLAS_OP_N,
                positions,
                static_cast<int>(out_channels),
                static_cast<int>(reduction),
                &alpha,
                columns.data_ptr<float>(), CUDA_R_32F, positions,
                weight.data_ptr<float>(), CUDA_R_32F,
                static_cast<int>(reduction),
                &beta,
                output_ptr, CUDA_R_32F, static_cast<int>(spatial),
                CUBLAS_COMPUTE_32F_PEDANTIC,
                CUBLAS_GEMM_DEFAULT);

            TORCH_CHECK(status == CUBLAS_STATUS_SUCCESS,
                        "cuBLAS GEMM failed: ", static_cast<int>(status));
        }
    }

    if (bias.numel() != 0) {
        TORCH_CHECK(bias.is_cuda() && bias.device() == x.device() &&
                    bias.scalar_type() == at::kFloat &&
                    bias.is_contiguous() && bias.numel() == out_channels,
                    "Invalid bias");

        const int blocks = static_cast<int>(
            std::min<int64_t>(4096, (output.numel() + 255) / 256));
        add_bias_kernel<<<blocks, 256, 0, stream>>>(
            output.data_ptr<float>(), bias.data_ptr<float>(),
            output.numel(), spatial, out_channels);
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    return output;
}
"""

_EXTENSION = None


def _get_extension():
    global _EXTENSION
    if _EXTENSION is None:
        _EXTENSION = load_inline(
            name="dilated_conv_unfold_cublas_v3",
            cpp_sources=_CPP_SOURCE,
            cuda_sources=_CUDA_SOURCE,
            functions=["dilated_conv_cuda"],
            extra_cflags=["-O3"],
            extra_cuda_cflags=["-O3"],
            extra_ldflags=["-lcublas"],
            verbose=False,
        )
    return _EXTENSION


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

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        conv = self.conv2d
        if (
            not x.is_cuda
            or x.ndim != 4
            or x.dtype != torch.float32
            or conv.weight.dtype != torch.float32
            or x.device != conv.weight.device
            or torch.is_autocast_enabled()
            or isinstance(conv.padding, str)
            or conv.weight[0].numel() > 65535
        ):
            return conv(x)

        extension = _get_extension()
        weight = conv.weight.contiguous()
        bias = (
            conv.bias.contiguous()
            if conv.bias is not None
            else x.new_empty(0)
        )

        return extension.dilated_conv_cuda(
            x.contiguous(),
            weight,
            bias,
            *conv.stride,
            *conv.padding,
            *conv.dilation,
        )
```