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

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

torch::Tensor linear_activation_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double multiplier,
    double negative_slope);
"""

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

__device__ __forceinline__ float activate(
    float value, float multiplier, float slope) {
    value *= multiplier;
    return value > 0.0f ? value : value * slope;
}

__global__ void activation_float4_kernel(
    float* data,
    int64_t size,
    float multiplier,
    float slope) {
    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    const int64_t stride =
        static_cast<int64_t>(blockDim.x) * gridDim.x;
    const int64_t vectors = size / 4;

    float4* vector_data = reinterpret_cast<float4*>(data);
    for (int64_t i = index; i < vectors; i += stride) {
        float4 value = vector_data[i];
        value.x = activate(value.x, multiplier, slope);
        value.y = activate(value.y, multiplier, slope);
        value.z = activate(value.z, multiplier, slope);
        value.w = activate(value.w, multiplier, slope);
        vector_data[i] = value;
    }

    const int64_t tail = vectors * 4 + index;
    if (tail < size) {
        data[tail] = activate(data[tail], multiplier, slope);
    }
}

template <typename scalar_t>
__global__ void activation_kernel(
    scalar_t* data,
    int64_t size,
    at::opmath_type<scalar_t> multiplier,
    at::opmath_type<scalar_t> slope) {
    using acc_t = at::opmath_type<scalar_t>;

    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    const int64_t stride =
        static_cast<int64_t>(blockDim.x) * gridDim.x;

    for (int64_t i = index; i < size; i += stride) {
        // Preserve the intermediate multiplication's dtype rounding.
        scalar_t multiplied =
            static_cast<scalar_t>(static_cast<acc_t>(data[i]) * multiplier);
        acc_t value = static_cast<acc_t>(multiplied);
        data[i] = static_cast<scalar_t>(
            value > acc_t(0) ? value : value * slope);
    }
}

torch::Tensor linear_activation_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double multiplier,
    double negative_slope) {
    const c10::cuda::CUDAGuard device_guard(x.device());

    // Retain PyTorch's optimized GEMM and fused bias addition.
    auto output = at::linear(x, weight, bias).contiguous();
    const int64_t size = output.numel();
    if (size == 0) {
        return output;
    }

    constexpr int threads = 256;
    const auto stream = at::cuda::getCurrentCUDAStream();

    if (output.scalar_type() == at::kFloat) {
        int64_t blocks = (size + threads * 4 - 1) / (threads * 4);
        blocks = blocks > 4096 ? 4096 : blocks;
        activation_float4_kernel<<<
            static_cast<unsigned int>(blocks), threads, 0, stream>>>(
                output.data_ptr<float>(),
                size,
                static_cast<float>(multiplier),
                static_cast<float>(negative_slope));
    } else {
        int64_t blocks = (size + threads - 1) / threads;
        blocks = blocks > 4096 ? 4096 : blocks;

        AT_DISPATCH_FLOATING_TYPES_AND2(
            at::ScalarType::Half,
            at::ScalarType::BFloat16,
            output.scalar_type(),
            "linear_activation_cuda",
            [&] {
                using acc_t = at::opmath_type<scalar_t>;
                activation_kernel<scalar_t><<<
                    static_cast<unsigned int>(blocks), threads, 0, stream>>>(
                        output.data_ptr<scalar_t>(),
                        size,
                        static_cast<acc_t>(multiplier),
                        static_cast<acc_t>(negative_slope));
            });
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = load_inline(
    name="linear_fused_scale_leaky_relu",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["linear_activation_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    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 not x.is_cuda:
            return self.leaky_relu(self.gemm(x) * self.multiplier)

        return _extension.linear_activation_cuda(
            x,
            self.gemm.weight,
            self.gemm.bias,
            float(self.multiplier),
            float(self.leaky_relu.negative_slope),
        )
```