```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 gemm_scale_leaky_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double multiplier,
    double negative_slope);
"""

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

__device__ __forceinline__ float activate(
    float x, float bias, float multiplier, float slope) {
    float y = __fmul_rn(__fadd_rn(x, bias), multiplier);
    return y > 0.0f ? y : __fmul_rn(y, slope);
}

__global__ void epilogue_vec4(
    float4* __restrict__ output,
    const float4* __restrict__ bias,
    int64_t count,
    int channels4,
    float multiplier,
    float slope) {
    int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
              + threadIdx.x;
    if (i >= count) return;

    float4 v = output[i];
    float4 b = bias[i % channels4];
    v.x = activate(v.x, b.x, multiplier, slope);
    v.y = activate(v.y, b.y, multiplier, slope);
    v.z = activate(v.z, b.z, multiplier, slope);
    v.w = activate(v.w, b.w, multiplier, slope);
    output[i] = v;
}

__global__ void epilogue_scalar(
    float* __restrict__ output,
    const float* __restrict__ bias,
    int64_t count,
    int channels,
    float multiplier,
    float slope) {
    int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
              + threadIdx.x;
    if (i < count) {
        output[i] = activate(
            output[i], bias[i % channels], multiplier, slope);
    }
}

torch::Tensor gemm_scale_leaky_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double multiplier,
    double negative_slope) {
    TORCH_CHECK(x.is_cuda() && weight.is_cuda() && bias.is_cuda(),
                "Expected CUDA tensors");
    TORCH_CHECK(x.device() == weight.device() &&
                x.device() == bias.device(),
                "All tensors must use the same device");
    TORCH_CHECK(x.scalar_type() == at::kFloat &&
                weight.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat,
                "Expected float32 tensors");
    TORCH_CHECK(x.dim() >= 1 && weight.dim() == 2 && bias.dim() == 1,
                "Invalid tensor dimensions");
    TORCH_CHECK(x.is_contiguous() && weight.is_contiguous() &&
                bias.is_contiguous(), "Expected contiguous tensors");

    const int64_t k64 = weight.size(1);
    const int64_t n64 = weight.size(0);
    TORCH_CHECK(k64 > 0 && x.size(-1) == k64 &&
                bias.numel() == n64, "Invalid feature dimensions");
    const int64_t m64 = x.numel() / k64;
    TORCH_CHECK(m64 <= INT_MAX && n64 <= INT_MAX && k64 <= INT_MAX,
                "GEMM dimensions exceed supported limits");

    const c10::cuda::CUDAGuard guard(x.device());
    auto sizes = x.sizes().vec();
    sizes.back() = n64;
    auto output = torch::empty(sizes, x.options());
    if (output.numel() == 0) return output;

    const int m = static_cast<int>(m64);
    const int n = static_cast<int>(n64);
    const int k = static_cast<int>(k64);

    auto stream = at::cuda::getCurrentCUDAStream();
    auto handle = at::cuda::getCurrentCUDABlasHandle();

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

    // Column-major interpretation computes output.T = weight @ x.T.
    // Explicit TF32 compute avoids depending on PyTorch's global setting.
    cublasStatus_t status = cublasGemmEx(
        handle,
        CUBLAS_OP_T,
        CUBLAS_OP_N,
        n, m, k,
        &alpha,
        weight.data_ptr<float>(), CUDA_R_32F, k,
        x.data_ptr<float>(), CUDA_R_32F, k,
        &beta,
        output.data_ptr<float>(), CUDA_R_32F, n,
        CUBLAS_COMPUTE_32F_FAST_TF32,
        CUBLAS_GEMM_DEFAULT_TENSOR_OP);

    TORCH_CHECK(status == CUBLAS_STATUS_SUCCESS,
                "cublasGemmEx failed: ", static_cast<int>(status));

    constexpr int threads = 256;
    const int64_t count = output.numel();
    if (n % 4 == 0 &&
        reinterpret_cast<uintptr_t>(bias.data_ptr<float>()) % 16 == 0) {
        const int64_t vectors = count / 4;
        epilogue_vec4<<<
            static_cast<unsigned int>((vectors + threads - 1) / threads),
            threads, 0, stream.stream()>>>(
                reinterpret_cast<float4*>(output.data_ptr<float>()),
                reinterpret_cast<const float4*>(bias.data_ptr<float>()),
                vectors, n / 4,
                static_cast<float>(multiplier),
                static_cast<float>(negative_slope));
    } else {
        epilogue_scalar<<<
            static_cast<unsigned int>((count + threads - 1) / threads),
            threads, 0, stream.stream()>>>(
                output.data_ptr<float>(), bias.data_ptr<float>(),
                count, n,
                static_cast<float>(multiplier),
                static_cast<float>(negative_slope));
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_ops = load_inline(
    name="tf32_gemm_inplace_scale_leaky_v2",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["gemm_scale_leaky_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    extra_ldflags=["-lcublas"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self, in_features, out_features, multiplier, negative_slope):
        super().__init__()
        self.gemm = nn.Linear(in_features, out_features)
        self.multiplier = multiplier
        self.leaky_relu = nn.LeakyReLU(negative_slope)

    def forward(self, x):
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and self.gemm.weight.dtype == torch.float32
            and self.gemm.bias.dtype == torch.float32
            and not torch.is_autocast_enabled()
            and self.gemm.in_features > 0
        ):
            return _ops.gemm_scale_leaky_cuda(
                x.contiguous(),
                self.gemm.weight.contiguous(),
                self.gemm.bias.contiguous(),
                float(self.multiplier),
                float(self.leaky_relu.negative_slope),
            )
        return self.leaky_relu(self.gemm(x) * self.multiplier)
```