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

cuda_src = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda_fp16.h>
#include <mma.h>

using namespace nvcuda;

#define BM 64
#define BN 128
#define BK 32
#define NTHREADS 256

// Rearranges ConvTranspose3d weight (IC, OC, KD, KH, KW) into the equivalent
// direct-convolution GEMM matrix A (OC, K) with K = IC*KD*KH*KW and flipped
// kernel indices, converted to fp16.
__global__ void prep_weight_kernel(const float* __restrict__ W,
                                   __half* __restrict__ A,
                                   int IC, int OC,
                                   int KD, int KH, int KW) {
    int KHW = KH * KW;
    int KKW = KW;
    int K = IC * KD * KHW;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= OC * K) return;
    int o = idx / K;
    int k = idx - o * K;
    int kw = k % KKW;
    int r = k / KKW;
    int kh = r % KH;
    r /= KH;
    int kd = r % KD;
    int i = r / KD;
    int src = (((i * OC + o) * KD + (KD - 1 - kd)) * KH + (KH - 1 - kh)) * KW + (KW - 1 - kw);
    A[idx] = __float2half(W[src]);
}

// Implicit GEMM ConvTranspose3d (stride=1, pad=0, out_pad=0, groups=1, dil=1).
// C[m, n] = sum_k A[m, k] * Xpad[k, n]
//   m = output channel, n = (batch, d, h, w) of output,
//   Xpad = input padded by (KD-1, KH-1, KW-1).
template <int KD, int KH, int KW>
__global__ void __launch_bounds__(NTHREADS)
ct3d_gemm_kernel(const __half* __restrict__ A,
                 const float* __restrict__ X,
                 float* __restrict__ Y,
                 const float* __restrict__ bias,
                 int Bt, int IC, int OC,
                 int ID, int IH, int IW,
                 int OD, int OH, int OW,
                 int K, int N) {
    constexpr int KHW = KH * KW;
    constexpr int KKW = KW;
    constexpr int KWK = KD * KH * KW;

    extern __shared__ char smem_raw[];
    __half* As0 = reinterpret_cast<__half*>(smem_raw);
    __half* Bs0 = As0 + BM * BK;
    __half* As1 = Bs0 + BK * BN;
    __half* Bs1 = As1 + BM * BK;
    float*  Cs  = reinterpret_cast<float*>(smem_raw);   // alias, used in epilogue

    const int tid  = threadIdx.x;
    const int bx   = blockIdx.x;
    const int by   = blockIdx.y;

    const int IS = ID * IH * IW;          // input channel stride
    const int OSP = OD * OH * OW;         // output spatial size

    // ---- per-thread n (column) assignment for B loads / epilogue ----
    const int n_loc = tid % BN;
    const int n_glob = bx * BN + n_loc;
    const bool valid_n = (n_glob < N);

    int x_base = 0;
    int od[KD], oh[KH], ow[KW];
    bool vd[KD], vh[KH], vw[KW];
    if (valid_n) {
        int w = n_glob % OW;
        int t1 = n_glob / OW;
        int h = t1 % OH;
        int t2 = t1 / OH;
        int d = t2 % OD;
        int b = t2 / OD;
        x_base = b * (IC * IS);
        #pragma unroll
        for (int kd = 0; kd < KD; ++kd) {
            int dd = d + kd - (KD - 1);
            vd[kd] = ((unsigned)dd < (unsigned)ID);
            od[kd] = dd * (IH * IW);
        }
        #pragma unroll
        for (int kh = 0; kh < KH; ++kh) {
            int hh = h + kh - (KH - 1);
            vh[kh] = ((unsigned)hh < (unsigned)IH);
            oh[kh] = hh * IW;
        }
        #pragma unroll
        for (int kw = 0; kw < KW; ++kw) {
            int ww = w + kw - (KW - 1);
            vw[kw] = ((unsigned)ww < (unsigned)IW);
            ow[kw] = ww;
        }
    }

    // A load assignment: thread t handles row m=t/4, 8 halves at k offset (t%4)*8
    const int a_m  = tid >> 2;
    const int a_kc = (tid & 3) * 8;
    // B load assignment: thread handles column n_loc, k rows [kg0, kg0+16)
    const int kg0 = (tid / BN) * 16;

    const int T = (K + BK - 1) / BK;

    // ---- lambda: load one K-tile into the given smem buffers ----
    auto load_tile = [&](int t, __half* As, __half* Bs) {
        int k0 = t * BK;
        // A tile
        if (a_m < OC) {
            if (k0 + a_kc + 8 <= K) {
                *reinterpret_cast<uint4*>(As + a_m * BK + a_kc) =
                    *reinterpret_cast<const uint4*>(A + (long)a_m * K + k0 + a_kc);
            } else {
                __half tmp[8];
                #pragma unroll
                for (int j = 0; j < 8; ++j) {
                    int kk = k0 + a_kc + j;
                    tmp[j] = (kk < K) ? A[(long)a_m * K + kk] : __float2half(0.f);
                }
                *reinterpret_cast<uint4*>(As + a_m * BK + a_kc) = *reinterpret_cast<uint4*>(tmp);
            }
        } else {
            *reinterpret_cast<uint4*>(As + a_m * BK + a_kc) = make_uint4(0, 0, 0, 0);
        }
        // B tile
        #pragma unroll 4
        for (int j = 0; j < 16; ++j) {
            int kg = k0 + kg0 + j;
            float v = 0.f;
            if (valid_n && kg < K) {
                int i  = kg / KWK;
                int r  = kg - i * KWK;
                int kd = r / KHW;
                int r2 = r - kd * KHW;
                int kh = r2 / KKW;
                int kw = r2 - kh * KKW;
                if (vd[kd] && vh[kh] && vw[kw]) {
                    v = X[x_base + i * IS + od[kd] + oh[kh] + ow[kw]];
                }
            }
            Bs[(kg0 + j) * BN + n_loc] = __float2half(v);
        }
    };

    // ---- accumulators ----
    wmma::fragment<wmma::accumulator, 16, 16, 16, float> c[2][2];
    #pragma unroll
    for (int i = 0; i < 2; ++i)
        #pragma unroll
        for (int j = 0; j < 2; ++j)
            wmma::fill_fragment(c[i][j], 0.0f);

    const int warp   = tid >> 5;
    const int warp_m = (warp >> 2) * 32;   // 2 warps along M
    const int warp_n = (warp & 3) * 32;    // 4 warps along N

    load_tile(0, As0, Bs0);
    __syncthreads();

    for (int t = 0; t < T; ++t) {
        __half* As = (t & 1) ? As1 : As0;
        __half* Bs = (t & 1) ? Bs1 : Bs0;

        // compute on current tile
        #pragma unroll
        for (int ks = 0; ks < BK / 16; ++ks) {
            wmma::fragment<wmma::matrix_a, 16, 16, 16, __half, wmma::row_major> af[2];
            wmma::fragment<wmma::matrix_b, 16, 16, 16, __half, wmma::row_major> bf[2];
            #pragma unroll
            for (int mi = 0; mi < 2; ++mi)
                wmma::load_matrix_sync(af[mi], As + (warp_m + mi * 16) * BK + ks * 16, BK);
            #pragma unroll
            for (int ni = 0; ni < 2; ++ni)
                wmma::load_matrix_sync(bf[ni], Bs + ks * 16 * BN + warp_n + ni * 16, BN);
            #pragma unroll
            for (int mi = 0; mi < 2; ++mi)
                #pragma unroll
                for (int ni = 0; ni < 2; ++ni)
                    wmma::mma_sync(c[mi][ni], af[mi], bf[ni], c[mi][ni]);
        }
        __syncthreads();
        // load next tile into the other buffer
        if (t + 1 < T) {
            __half* AsN = (t & 1) ? As0 : As1;
            __half* BsN = (t & 1) ? Bs0 : Bs1;
            load_tile(t + 1, AsN, BsN);
            __syncthreads();
        }
    }

    // ---- epilogue: stage through smem ----
    __syncthreads();
    #pragma unroll
    for (int mi = 0; mi < 2; ++mi)
        #pragma unroll
        for (int ni = 0; ni < 2; ++ni)
            wmma::store_matrix_sync(Cs + (warp_m + mi * 16) * BN + warp_n + ni * 16,
                                    c[mi][ni], BN, wmma::mem_row_major);
    __syncthreads();

    // write out: thread's n stays fixed; m advances by 2 per iteration
    {
        int b = 0, sp = 0;
        if (valid_n) { b = n_glob / OSP; sp = n_glob - b * OSP; }
        int m0 = tid / BN;                  // 0 or 1
        int o = by * BM + m0;
        int addr = (b * OC + o) * OSP + sp;
        #pragma unroll 4
        for (int r = 0; r < BM / 2; ++r) {
            if (valid_n && o < OC) {
                float v = Cs[(m0 + 2 * r) * BN + n_loc];
                if (bias != nullptr) v += bias[o];
                Y[addr] = v;
            }
            o += 2;
            addr += 2 * OSP;
        }
    }
}

torch::Tensor conv_transpose3d_forward(torch::Tensor x, torch::Tensor weight,
                                       torch::Tensor bias,
                                       int64_t KD, int64_t KH, int64_t KW) {
    TORCH_CHECK(x.is_cuda() && x.dtype() == torch::kFloat);
    TORCH_CHECK(x.dim() == 5 && weight.dim() == 5);

    auto xc = x.contiguous();
    auto wc = weight.contiguous();

    const int Bt = xc.size(0);
    const int IC = xc.size(1);
    const int ID = xc.size(2);
    const int IH = xc.size(3);
    const int IW = xc.size(4);
    const int OC = wc.size(1);   // weight: (IC, OC, KD, KH, KW)

    const int OD = ID + (int)KD - 1;
    const int OH = IH + (int)KH - 1;
    const int OW = IW + (int)KW - 1;

    auto Y = torch::empty({Bt, OC, OD, OH, OW}, xc.options());

    const int K = IC * (int)(KD * KH * KW);
    const long N = (long)Bt * OD * OH * OW;

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

    // weight rearrange -> A (OC, K) in fp16
    auto Ahalf = torch::empty({(long)OC * K}, xc.options().dtype(torch::kHalf));
    {
        int total = OC * K;
        int threads = 256;
        int blocks = (total + threads - 1) / threads;
        prep_weight_kernel<<<blocks, threads, 0, stream>>>(
            wc.data_ptr<float>(),
            reinterpret_cast<__half*>(Ahalf.data_ptr<at::Half>()),
            IC, OC, (int)KD, (int)KH, (int)KW);
    }

    dim3 grid((unsigned)((N + BN - 1) / BN), (unsigned)((OC + BM - 1) / BM));
    const int smem_bytes = BM * BN * (int)sizeof(float);   // 32 KB (covers double-buffered A/B too)

    const float* bias_ptr = nullptr;
    torch::Tensor biasc;
    if (bias.defined() && bias.numel() > 0) {
        biasc = bias.contiguous();
        bias_ptr = biasc.data_ptr<float>();
    }

    auto launch = [&](auto kd, auto kh, auto kw) {
        constexpr int KDc = decltype(kd)::value;
        constexpr int KHc = decltype(kh)::value;
        constexpr int KWc = decltype(kw)::value;
        auto kfn = ct3d_gemm_kernel<KDc, KHc, KWc>;
        static bool configured = false;
        if (!configured) {
            cudaFuncSetAttribute(kfn, cudaFuncAttributeMaxDynamicSharedMemorySize, smem_bytes);
            configured = true;
        }
        kfn<<<grid, NTHREADS, smem_bytes, stream>>>(
            reinterpret_cast<const __half*>(Ahalf.data_ptr<at::Half>()),
            xc.data_ptr<float>(), Y.data_ptr<float>(), bias_ptr,
            Bt, IC, OC, ID, IH, IW, OD, OH, OW, K, (int)N);
    };

    if (KD == 3 && KH == 5 && KW == 5) {
        launch(std::integral_constant<int, 3>{}, std::integral_constant<int, 5>{}, std::integral_constant<int, 5>{});
    } else if (KD == 3 && KH == 3 && KW == 3) {
        launch(std::integral_constant<int, 3>{}, std::integral_constant<int, 3>{}, std::integral_constant<int, 3>{});
    } else if (KD == 1 && KH == 1 && KW == 1) {
        launch(std::integral_constant<int, 1>{}, std::integral_constant<int, 1>{}, std::integral_constant<int, 1>{});
    } else if (KD == 5 && KH == 5 && KW == 5) {
        launch(std::integral_constant<int, 5>{}, std::integral_constant<int, 5>{}, std::integral_constant<int, 5>{});
    } else if (KD == 2 && KH == 2 && KW == 2) {
        launch(std::integral_constant<int, 2>{}, std::integral_constant<int, 2>{}, std::integral_constant<int, 2>{});
    } else if (KD == 4 && KH == 4 && KW == 4) {
        launch(std::integral_constant<int, 4>{}, std::integral_constant<int, 4>{}, std::integral_constant<int, 4>{});
    } else if (KD == 7 && KH == 7 && KW == 7) {
        launch(std::integral_constant<int, 7>{}, std::integral_constant<int, 7>{}, std::integral_constant<int, 7>{});
    } else if (KD == 1 && KH == 3 && KW == 3) {
        launch(std::integral_constant<int, 1>{}, std::integral_constant<int, 3>{}, std::integral_constant<int, 3>{});
    } else if (KD == 3 && KH == 1 && KW == 1) {
        launch(std::integral_constant<int, 3>{}, std::integral_constant<int, 1>{}, std::integral_constant<int, 1>{});
    } else if (KD == 1 && KH == 5 && KW == 5) {
        launch(std::integral_constant<int, 1>{}, std::integral_constant<int, 5>{}, std::integral_constant<int, 5>{});
    } else if (KD == 5 && KH == 3 && KW == 3) {
        launch(std::integral_constant<int, 5>{}, std::integral_constant<int, 3>{}, std::integral_constant<int, 3>{});
    } else if (KD == 3 && KH == 7 && KW == 7) {
        launch(std::integral_constant<int, 3>{}, std::integral_constant<int, 7>{}, std::integral_constant<int, 7>{});
    } else {
        TORCH_CHECK(false, "Unsupported kernel size for custom ConvTranspose3d fast path");
    }

    return Y;
}
"""

cpp_src = r"""
torch::Tensor conv_transpose3d_forward(torch::Tensor x, torch::Tensor weight,
                                       torch::Tensor bias,
                                       int64_t KD, int64_t KH, int64_t KW);
"""

_mod = load_inline(
    name="conv_transpose3d_fast",
    cpp_sources=cpp_src,
    cuda_sources=cuda_src,
    functions=["conv_transpose3d_forward"],
    verbose=False,
    extra_cuda_cflags=["-O3", "--use_fast_math", "-std=c++17"],
)


class ModelNew(nn.Module):
    """
    ConvTranspose3d accelerated with a custom implicit-GEMM CUDA kernel that
    uses fp16 tensor cores with fp32 accumulation. Falls back to the native
    PyTorch implementation for configurations outside the fast path.
    """

    def __init__(self, in_channels: int, out_channels: int, kernel_size: tuple,
                 stride: tuple = (1, 1, 1), padding: tuple = (0, 0, 0),
                 output_padding: tuple = (0, 0, 0), groups: int = 1,
                 bias: bool = False):
        super(ModelNew, self).__init__()
        self.conv_transpose3d = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size, stride=stride,
            padding=padding, output_padding=output_padding, groups=groups,
            bias=bias)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        conv = self.conv_transpose3d
        if (x.is_cuda and x.dtype == torch.float32
                and tuple(conv.stride) == (1, 1, 1)
                and tuple(conv.padding) == (0, 0, 0)
                and tuple(conv.output_padding) == (0, 0, 0)
                and tuple(conv.dilation) == (1, 1, 1)
                and conv.groups == 1):
            kd, kh, kw = conv.kernel_size
            bias_t = conv.bias if conv.bias is not None else torch.empty(
                0, device=x.device, dtype=x.dtype)
            return _mod.conv_transpose3d_forward(
                x, conv.weight, bias_t, kd, kh, kw)
        return self.conv_transpose3d(x)
```