```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 depthwise_conv_cuda(
    torch::Tensor x, torch::Tensor weight, torch::Tensor bias,
    bool has_bias, int64_t stride_h, int64_t stride_w,
    int64_t padding_h, int64_t padding_w,
    int64_t dilation_h, int64_t dilation_w);
"""

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

__global__ void depthwise_3x7_float(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int64_t channels,
    int64_t height,
    int64_t width,
    int64_t out_height,
    int64_t out_width)
{
    constexpr int ROWS = 4;
    const int64_t ox = blockIdx.x * blockDim.x + threadIdx.x;
    const int64_t oy = blockIdx.y * ROWS;
    const int64_t plane = blockIdx.z;
    if (ox >= out_width) return;

    const int64_t channel = plane % channels;
    const float* input = x + plane * height * width;
    const float* filter = weight + channel * 21;

    float w[21];
    #pragma unroll
    for (int k = 0; k < 21; ++k)
        w[k] = filter[k];

    float accum[ROWS];
    const float initial = bias ? bias[channel] : 0.0f;
    #pragma unroll
    for (int q = 0; q < ROWS; ++q)
        accum[q] = initial;

    // Reuse each input value across up to three neighboring output rows.
    #pragma unroll
    for (int r = 0; r < ROWS + 2; ++r) {
        #pragma unroll
        for (int k = 0; k < 7; ++k) {
            const float value = oy + r < height
                ? input[(oy + r) * width + ox + k]
                : 0.0f;

            #pragma unroll
            for (int q = 0; q < ROWS; ++q) {
                if (r >= q && r < q + 3)
                    accum[q] = fmaf(value, w[(r - q) * 7 + k], accum[q]);
            }
        }
    }

    float* destination = output + plane * out_height * out_width;
    #pragma unroll
    for (int q = 0; q < ROWS; ++q) {
        if (oy + q < out_height)
            destination[(oy + q) * out_width + ox] = accum[q];
    }
}

template <typename scalar_t>
__global__ void depthwise_generic(
    const scalar_t* __restrict__ x,
    const scalar_t* __restrict__ weight,
    const scalar_t* __restrict__ bias,
    scalar_t* __restrict__ output,
    int64_t total,
    int64_t channels,
    int64_t height,
    int64_t width,
    int64_t out_height,
    int64_t out_width,
    int64_t kernel_h,
    int64_t kernel_w,
    int64_t stride_h,
    int64_t stride_w,
    int64_t padding_h,
    int64_t padding_w,
    int64_t dilation_h,
    int64_t dilation_w)
{
    using acc_t = at::acc_type<scalar_t, true>;

    for (int64_t index = static_cast<int64_t>(blockIdx.x) * blockDim.x
                          + threadIdx.x;
         index < total;
         index += static_cast<int64_t>(blockDim.x) * gridDim.x) {
        const int64_t ox = index % out_width;
        const int64_t oy = (index / out_width) % out_height;
        const int64_t plane = index / (out_height * out_width);
        const int64_t channel = plane % channels;
        const int64_t iy0 = oy * stride_h - padding_h;
        const int64_t ix0 = ox * stride_w - padding_w;

        const scalar_t* input = x + plane * height * width;
        const scalar_t* filter = weight + channel * kernel_h * kernel_w;
        acc_t sum = acc_t(0);

        for (int64_t kh = 0; kh < kernel_h; ++kh) {
            const int64_t iy = iy0 + kh * dilation_h;
            if (iy < 0 || iy >= height) continue;

            for (int64_t kw = 0; kw < kernel_w; ++kw) {
                const int64_t ix = ix0 + kw * dilation_w;
                if (ix >= 0 && ix < width) {
                    sum += static_cast<acc_t>(input[iy * width + ix])
                         * static_cast<acc_t>(filter[kh * kernel_w + kw]);
                }
            }
        }

        if (bias) sum += static_cast<acc_t>(bias[channel]);
        output[index] = static_cast<scalar_t>(sum);
    }
}

torch::Tensor depthwise_conv_cuda(
    torch::Tensor x, torch::Tensor weight, torch::Tensor bias,
    bool has_bias, int64_t stride_h, int64_t stride_w,
    int64_t padding_h, int64_t padding_w,
    int64_t dilation_h, int64_t dilation_w)
{
    TORCH_CHECK(x.is_cuda() && weight.is_cuda(), "Expected CUDA tensors");
    TORCH_CHECK(x.dim() == 4, "Expected NCHW input");
    TORCH_CHECK(weight.dim() == 4 && weight.size(1) == 1,
                "Expected depthwise convolution weights");
    TORCH_CHECK(x.device() == weight.device(), "Device mismatch");
    TORCH_CHECK(x.scalar_type() == weight.scalar_type(), "Dtype mismatch");
    TORCH_CHECK(x.size(1) == weight.size(0), "Channel mismatch");
    TORCH_CHECK(stride_h > 0 && stride_w > 0 &&
                dilation_h > 0 && dilation_w > 0 &&
                padding_h >= 0 && padding_w >= 0,
                "Invalid convolution parameters");

    const c10::cuda::CUDAGuard guard(x.device());
    x = x.contiguous();
    weight = weight.contiguous();

    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 kernel_h = weight.size(2);
    const int64_t kernel_w = weight.size(3);

    if (has_bias) {
        TORCH_CHECK(bias.device() == x.device() &&
                    bias.scalar_type() == x.scalar_type() &&
                    bias.numel() == channels, "Invalid bias");
        bias = bias.contiguous();
    }

    const int64_t available_h =
        height + 2 * padding_h - dilation_h * (kernel_h - 1) - 1;
    const int64_t available_w =
        width + 2 * padding_w - dilation_w * (kernel_w - 1) - 1;
    TORCH_CHECK(available_h >= 0 && available_w >= 0,
                "Kernel is larger than the padded input");

    const int64_t out_height = available_h / stride_h + 1;
    const int64_t out_width = available_w / stride_w + 1;
    auto output = torch::empty(
        {batch, channels, out_height, out_width}, x.options());

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

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (x.scalar_type() == at::kFloat &&
        kernel_h == 3 && kernel_w == 7 &&
        stride_h == 1 && stride_w == 1 &&
        padding_h == 0 && padding_w == 0 &&
        dilation_h == 1 && dilation_w == 1 &&
        batch * channels <= 65535 &&
        (out_height + 3) / 4 <= 65535) {
        const dim3 block(256);
        const dim3 grid(
            static_cast<unsigned int>((out_width + 255) / 256),
            static_cast<unsigned int>((out_height + 3) / 4),
            static_cast<unsigned int>(batch * channels));

        depthwise_3x7_float<<<grid, block, 0, stream>>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            has_bias ? bias.data_ptr<float>() : nullptr,
            output.data_ptr<float>(),
            channels, height, width, out_height, out_width);
    } else {
        const int64_t total = output.numel();
        const int blocks = static_cast<int>(
            std::min<int64_t>((total + 255) / 256, 65535));

        AT_DISPATCH_FLOATING_TYPES_AND2(
            at::ScalarType::Half, at::ScalarType::BFloat16,
            x.scalar_type(), "depthwise_generic_cuda", [&] {
                depthwise_generic<scalar_t><<<blocks, 256, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    weight.data_ptr<scalar_t>(),
                    has_bias ? bias.data_ptr<scalar_t>() : nullptr,
                    output.data_ptr<scalar_t>(),
                    total, channels, height, width, out_height, out_width,
                    kernel_h, kernel_w, stride_h, stride_w,
                    padding_h, padding_w, dilation_h, dilation_w);
            });
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_depthwise_cuda = load_inline(
    name="depthwise_asymmetric_row_reuse_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["depthwise_conv_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_h: int,
        kernel_size_w: int,
        stride_h: int = 1,
        stride_w: int = 1,
        padding_h: int = 0,
        padding_w: int = 0,
        dilation_h: int = 1,
        dilation_w: int = 1,
        groups: int = 1,
        bias: bool = False,
    ):
        super().__init__()
        self.conv2d = nn.Conv2d(
            in_channels,
            in_channels,
            (kernel_size_h, kernel_size_w),
            stride=(stride_h, stride_w),
            padding=(padding_h, padding_w),
            dilation=(dilation_h, dilation_w),
            groups=in_channels,
            bias=bias,
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        if not x.is_cuda:
            return self.conv2d(x)

        conv = self.conv2d
        return _depthwise_cuda.depthwise_conv_cuda(
            x,
            conv.weight,
            conv.bias if conv.bias is not None else conv.weight,
            conv.bias is not None,
            conv.stride[0],
            conv.stride[1],
            conv.padding[0],
            conv.padding[1],
            conv.dilation[0],
            conv.dilation[1],
        )
```