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

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

__global__ void depthwise3x3_kernel(
    const float* __restrict__ input,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int channels,
    int height,
    int width,
    int out_height,
    int out_width,
    int tiles_x,
    int tiles_y)
{
    const int tile_x = blockIdx.x % tiles_x;
    const int remainder = blockIdx.x / tiles_x;
    const int tile_y = remainder % tiles_y;
    const int plane = remainder / tiles_y;
    const int channel = plane % channels;

    const int x = tile_x * 512 + threadIdx.x * 4;
    const int y = tile_y * 4;
    if (x >= out_width) return;

    const float* src =
        input + static_cast<int64_t>(plane) * height * width;
    float* dst =
        output + static_cast<int64_t>(plane) * out_height * out_width;

    float k[9];
    #pragma unroll
    for (int i = 0; i < 9; ++i)
        k[i] = weight[channel * 9 + i];

    float acc[4][4];
    #pragma unroll
    for (int oy = 0; oy < 4; ++oy) {
        #pragma unroll
        for (int ox = 0; ox < 4; ++ox)
            acc[oy][ox] = 0.0f;
    }

    // Reuse each input row across up to three output rows.
    #pragma unroll
    for (int row = 0; row < 6; ++row) {
        if (y + row < height) {
            const float* p =
                src + static_cast<int64_t>(y + row) * width + x;
            const float4 first = *reinterpret_cast<const float4*>(p);

            float2 last = make_float2(0.0f, 0.0f);
            if (x + 4 < width)
                last = *reinterpret_cast<const float2*>(p + 4);

            float v[6] = {
                first.x, first.y, first.z, first.w, last.x, last.y
            };

            #pragma unroll
            for (int oy = 0; oy < 4; ++oy) {
                if (row >= oy && row < oy + 3) {
                    const int kr = (row - oy) * 3;
                    #pragma unroll
                    for (int ox = 0; ox < 4; ++ox) {
                        acc[oy][ox] = fmaf(v[ox],     k[kr],     acc[oy][ox]);
                        acc[oy][ox] = fmaf(v[ox + 1], k[kr + 1], acc[oy][ox]);
                        acc[oy][ox] = fmaf(v[ox + 2], k[kr + 2], acc[oy][ox]);
                    }
                }
            }
        }
    }

    const float b = bias ? bias[channel] : 0.0f;
    #pragma unroll
    for (int oy = 0; oy < 4; ++oy) {
        if (y + oy < out_height) {
            float* p = dst + static_cast<int64_t>(y + oy) * out_width + x;
            *reinterpret_cast<float2*>(p) =
                make_float2(acc[oy][0] + b, acc[oy][1] + b);
            if (x + 2 < out_width) {
                *reinterpret_cast<float2*>(p + 2) =
                    make_float2(acc[oy][2] + b, acc[oy][3] + b);
            }
        }
    }
}

torch::Tensor depthwise3x3_cuda(
    torch::Tensor input,
    torch::Tensor weight,
    c10::optional<torch::Tensor> bias)
{
    const c10::cuda::CUDAGuard device_guard(input.device());

    const int batch = input.size(0);
    const int channels = input.size(1);
    const int height = input.size(2);
    const int width = input.size(3);
    const int out_height = height - 2;
    const int out_width = width - 2;

    auto output = torch::empty(
        {batch, channels, out_height, out_width}, input.options());
    if (output.numel() == 0)
        return output;

    const int tiles_x = (out_width + 511) / 512;
    const int tiles_y = (out_height + 3) / 4;
    const int64_t blocks =
        static_cast<int64_t>(batch) * channels * tiles_x * tiles_y;

    TORCH_CHECK(blocks <= 2147483647LL, "Input exceeds supported grid size");

    const float* bias_ptr =
        bias.has_value() ? bias->data_ptr<float>() : nullptr;

    depthwise3x3_kernel<<<
        static_cast<unsigned int>(blocks), 128, 0,
        at::cuda::getCurrentCUDAStream()
    >>>(
        input.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias_ptr,
        output.data_ptr<float>(),
        channels, height, width, out_height, out_width, tiles_x, tiles_y
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

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

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

_depthwise_ext = load_inline(
    name="depthwise_3x3_tiled_cuda_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["depthwise3x3_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


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

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        conv = self.conv2d
        weight = conv.weight
        bias = conv.bias

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 4
            and x.is_contiguous()
            and x.shape[1] == conv.in_channels
            and x.shape[2] >= 3
            and x.shape[3] >= 4
            and x.shape[3] % 4 == 0
            and conv.in_channels == conv.out_channels
            and conv.kernel_size == (3, 3)
            and conv.stride == (1, 1)
            and conv.padding == (0, 0)
            and conv.dilation == (1, 1)
            and weight.dtype == torch.float32
            and weight.device == x.device
            and weight.is_contiguous()
            and (
                bias is None
                or (
                    bias.dtype == torch.float32
                    and bias.device == x.device
                    and bias.is_contiguous()
                )
            )
        ):
            return _depthwise_ext.depthwise3x3_cuda(x, weight, bias)

        return conv(x)
```