```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>
#include <algorithm>
#include <cstdint>

__device__ __forceinline__ float pool_max(float a, float b) {
    return (b > a || isnan(b)) ? b : a;
}

template <bool Interior>
__device__ __forceinline__ float spatial_max(
    const float* __restrict__ plane,
    int height,
    int width,
    int y,
    int x
) {
    float result = -CUDART_INF_F;

    #pragma unroll
    for (int ky = 0; ky < 3; ++ky) {
        const int iy = y + 3 * ky;

        #pragma unroll
        for (int kx = 0; kx < 3; ++kx) {
            const int ix = x + 3 * kx;
            if (Interior ||
                (iy >= 0 && iy < height && ix >= 0 && ix < width)) {
                result = pool_max(
                    result, plane[static_cast<int64_t>(iy) * width + ix]
                );
            }
        }
    }
    return result;
}

__global__ void maxpool3d_3_2_1_3_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    int depth,
    int height,
    int width,
    int out_depth,
    int out_height,
    int out_width,
    int64_t channel_count,
    int pairs_per_channel
) {
    const int ox = blockIdx.x * blockDim.x + threadIdx.x;
    const int oy = blockIdx.y * blockDim.y + threadIdx.y;

    if (ox >= out_width || oy >= out_height)
        return;

    const int ix = 2 * ox - 1;
    const int iy = 2 * oy - 1;
    const bool interior =
        ix >= 0 && ix + 6 < width &&
        iy >= 0 && iy + 6 < height;

    const int64_t input_plane = static_cast<int64_t>(height) * width;
    const int64_t output_plane =
        static_cast<int64_t>(out_height) * out_width;
    const int64_t total_pairs = channel_count * pairs_per_channel;

    for (int64_t pair = blockIdx.z;
         pair < total_pairs;
         pair += gridDim.z) {
        const int64_t channel = pair / pairs_per_channel;
        const int local_pair = pair % pairs_per_channel;

        // Pair output depths separated by three: their windows share a plane.
        const int oz = (local_pair / 3) * 6 + local_pair % 3;
        if (oz >= out_depth)
            continue;

        const bool second_valid = oz + 3 < out_depth;
        const int iz = 2 * oz - 1;
        const float* channel_input =
            input + channel * depth * input_plane;

        float first = -CUDART_INF_F;
        float second = -CUDART_INF_F;

        #pragma unroll
        for (int plane_id = 0; plane_id < 5; ++plane_id) {
            if (plane_id < 3 || second_valid) {
                const int z = iz + 3 * plane_id;
                float value = -CUDART_INF_F;

                if (z >= 0 && z < depth) {
                    const float* plane =
                        channel_input + static_cast<int64_t>(z) * input_plane;
                    if (interior) {
                        value = spatial_max<true>(
                            plane, height, width, iy, ix
                        );
                    } else {
                        value = spatial_max<false>(
                            plane, height, width, iy, ix
                        );
                    }
                }

                if (plane_id <= 2)
                    first = pool_max(first, value);
                if (plane_id >= 2)
                    second = pool_max(second, value);
            }
        }

        const int64_t offset =
            (channel * out_depth + oz) * output_plane +
            static_cast<int64_t>(oy) * out_width + ox;

        output[offset] = first;
        if (second_valid)
            output[offset + 3 * output_plane] = second;
    }
}

torch::Tensor maxpool3d_specialized(torch::Tensor input) {
    TORCH_CHECK(input.is_cuda(), "Expected a CUDA tensor");
    TORCH_CHECK(input.scalar_type() == at::kFloat, "Expected float32");
    TORCH_CHECK(input.is_contiguous(), "Expected contiguous input");
    TORCH_CHECK(input.dim() == 5, "Expected a five-dimensional input");

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

    const int depth = input.size(2);
    const int height = input.size(3);
    const int width = input.size(4);

    TORCH_CHECK(
        depth >= 5 && height >= 5 && width >= 5,
        "Input spatial dimensions must be at least five"
    );

    const int out_depth = (depth - 5) / 2 + 1;
    const int out_height = (height - 5) / 2 + 1;
    const int out_width = (width - 5) / 2 + 1;

    auto output = torch::empty(
        {input.size(0), input.size(1), out_depth, out_height, out_width},
        input.options()
    );

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

    const int64_t channel_count = input.size(0) * input.size(1);
    const int pairs_per_channel = ((out_depth + 5) / 6) * 3;

    const dim3 block(32, 8);
    const dim3 grid(
        (out_width + 31) / 32,
        (out_height + 7) / 8,
        static_cast<unsigned int>(std::min<int64_t>(
            channel_count * pairs_per_channel, 65535
        ))
    );

    maxpool3d_3_2_1_3_kernel<<<
        grid, block, 0, at::cuda::getCurrentCUDAStream()
    >>>(
        input.data_ptr<float>(),
        output.data_ptr<float>(),
        depth,
        height,
        width,
        out_depth,
        out_height,
        out_width,
        channel_count,
        pairs_per_channel
    );

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = load_inline(
    name="maxpool3d_paired_depth_cuda",
    cpp_sources=(
        "torch::Tensor maxpool3d_specialized(torch::Tensor input);"
    ),
    cuda_sources=_cuda_source,
    functions=["maxpool3d_specialized"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


def _triple(value):
    if isinstance(value, int):
        return (value, value, value)
    return tuple(value)


class ModelNew(nn.Module):
    def __init__(
        self,
        kernel_size: int,
        stride: int = None,
        padding: int = 0,
        dilation: int = 1,
        return_indices: bool = False,
        ceil_mode: bool = False,
    ):
        super().__init__()
        self.maxpool = nn.MaxPool3d(
            kernel_size=kernel_size,
            stride=stride,
            padding=padding,
            dilation=dilation,
            return_indices=return_indices,
            ceil_mode=ceil_mode,
        )
        effective_stride = kernel_size if stride is None else stride
        self._specialized = (
            _triple(kernel_size) == (3, 3, 3)
            and _triple(effective_stride) == (2, 2, 2)
            and _triple(padding) == (1, 1, 1)
            and _triple(dilation) == (3, 3, 3)
            and not return_indices
            and not ceil_mode
        )

    def forward(self, x: torch.Tensor):
        if (
            self._specialized
            and x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 5
            and x.is_contiguous()
            and min(x.shape[2:]) >= 5
            and x.shape[1] > 0
            and (x.shape[3] - 3) // 2 <= 65535 * 8
            and not (torch.is_grad_enabled() and x.requires_grad)
        ):
            return _extension.maxpool3d_specialized(x)
        return self.maxpool(x)
```