```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 fused_epilogue_cuda(
    torch::Tensor x, double subtract_value, double multiply_value);
"""

_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>
#include <algorithm>

__device__ __forceinline__ float transform(
    float x, float subtract_value, float multiply_value) {
    float y = __fmul_rn(__fsub_rn(x, subtract_value), multiply_value);
    return y <= 0.0f ? 0.0f : y;
}

__global__ void epilogue_float4(
    float4* x,
    int64_t count,
    float subtract_value,
    float multiply_value) {
    const int64_t stride = static_cast<int64_t>(blockDim.x) * gridDim.x;
    for (int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
                     + threadIdx.x;
         i < count; i += stride) {
        float4 v = x[i];
        v.x = transform(v.x, subtract_value, multiply_value);
        v.y = transform(v.y, subtract_value, multiply_value);
        v.z = transform(v.z, subtract_value, multiply_value);
        v.w = transform(v.w, subtract_value, multiply_value);
        x[i] = v;
    }
}

template <typename scalar_t>
__global__ void epilogue_scalar(
    scalar_t* x,
    int64_t count,
    at::acc_type<scalar_t, true> subtract_value,
    at::acc_type<scalar_t, true> multiply_value) {
    using acc_t = at::acc_type<scalar_t, true>;
    const int64_t stride = static_cast<int64_t>(blockDim.x) * gridDim.x;

    for (int64_t i = static_cast<int64_t>(blockIdx.x) * blockDim.x
                     + threadIdx.x;
         i < count; i += stride) {
        // Preserve intermediate rounding for half and bfloat16 tensors.
        scalar_t subtracted =
            static_cast<scalar_t>(static_cast<acc_t>(x[i]) - subtract_value);
        scalar_t multiplied =
            static_cast<scalar_t>(static_cast<acc_t>(subtracted) * multiply_value);
        x[i] = static_cast<acc_t>(multiplied) <= acc_t(0)
            ? static_cast<scalar_t>(0)
            : multiplied;
    }
}

torch::Tensor fused_epilogue_cuda(
    torch::Tensor x, double subtract_value, double multiply_value) {
    TORCH_CHECK(x.is_cuda(), "Expected a CUDA tensor");
    TORCH_CHECK(x.is_contiguous(), "Expected a contiguous tensor");

    const c10::cuda::CUDAGuard device_guard(x.device());
    const int64_t count = x.numel();
    if (count == 0) {
        return x;
    }

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

    if (x.scalar_type() == at::kFloat && count % 4 == 0) {
        const int64_t vector_count = count / 4;
        const int blocks = static_cast<int>(
            std::min<int64_t>((vector_count + threads - 1) / threads, 65535));

        epilogue_float4<<<blocks, threads, 0, stream>>>(
            reinterpret_cast<float4*>(x.data_ptr<float>()),
            vector_count,
            static_cast<float>(subtract_value),
            static_cast<float>(multiply_value));
    } else {
        const int blocks = static_cast<int>(
            std::min<int64_t>((count + threads - 1) / threads, 65535));

        AT_DISPATCH_FLOATING_TYPES_AND2(
            at::ScalarType::Half,
            at::ScalarType::BFloat16,
            x.scalar_type(),
            "fused_epilogue_cuda",
            [&] {
                using acc_t = at::acc_type<scalar_t, true>;
                epilogue_scalar<scalar_t><<<blocks, threads, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    count,
                    static_cast<acc_t>(subtract_value),
                    static_cast<acc_t>(multiply_value));
            });
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return x;
}
"""

_fused_ops = load_inline(
    name="linear_sub_mul_relu_epilogue",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_epilogue_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3", "--fmad=false"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_features, out_features, subtract_value, multiply_value
    ):
        super().__init__()
        self.linear = nn.Linear(in_features, out_features)
        self.subtract_value = subtract_value
        self.multiply_value = multiply_value

    def forward(self, x):
        y = self.linear(x)

        # Retain autograd support when gradients are requested.
        if not y.is_cuda or (torch.is_grad_enabled() and y.requires_grad):
            return torch.relu(
                (y - self.subtract_value) * self.multiply_value
            )

        return _fused_ops.fused_epilogue_cuda(
            y.contiguous(), self.subtract_value, self.multiply_value
        )
```