```python
import math
import torch
import torch.nn as nn
import torch.nn.functional as F
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 <cuda_fp16.h>
#include <cublasLt.h>
#include <algorithm>
#include <cstdint>
#include <map>
#include <memory>
#include <tuple>

static void check_lt(cublasStatus_t status) {
    TORCH_CHECK(status == CUBLAS_STATUS_SUCCESS,
                "cuBLASLt error: ", static_cast<int>(status));
}

__global__ void pack_weight_kernel(
    const float* __restrict__ src,
    half* __restrict__ dst,
    int64_t vectors
) {
    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < vectors;
         i += int64_t(blockDim.x) * gridDim.x) {
        float4 v = reinterpret_cast<const float4*>(src)[i];
        half2* d = reinterpret_cast<half2*>(dst) + 2 * i;
        d[0] = __floats2half2_rn(v.x * 256.0f, v.y * 256.0f);
        d[1] = __floats2half2_rn(v.z * 256.0f, v.w * 256.0f);
    }
}

__global__ void pack_input_bias_kernel(
    const float* __restrict__ src,
    half* __restrict__ dst,
    const float* __restrict__ bias,
    float* __restrict__ adjusted_bias,
    int64_t vectors,
    int bias_vectors,
    float sub,
    float mul
) {
    const int64_t first =
        int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    const int64_t step = int64_t(blockDim.x) * gridDim.x;

    for (int64_t i = first; i < vectors; i += step) {
        float4 v = reinterpret_cast<const float4*>(src)[i];
        half2* d = reinterpret_cast<half2*>(dst) + 2 * i;
        d[0] = __floats2half2_rn(v.x, v.y);
        d[1] = __floats2half2_rn(v.z, v.w);
    }

    for (int64_t i = first; i < bias_vectors; i += step) {
        float4 b = reinterpret_cast<const float4*>(bias)[i];
        b.x = __fmul_rn(__fsub_rn(b.x, sub), mul);
        b.y = __fmul_rn(__fsub_rn(b.y, sub), mul);
        b.z = __fmul_rn(__fsub_rn(b.z, sub), mul);
        b.w = __fmul_rn(__fsub_rn(b.w, sub), mul);
        reinterpret_cast<float4*>(adjusted_bias)[i] = b;
    }
}

torch::Tensor pack_weight_cuda(torch::Tensor weight) {
    const c10::cuda::CUDAGuard guard(weight.device());
    weight = weight.contiguous();

    auto packed = torch::empty(
        weight.sizes(), weight.options().dtype(at::kHalf));
    const int64_t vectors = weight.numel() / 4;
    const int blocks = static_cast<int>(
        std::min<int64_t>((vectors + 255) / 256, 8192));

    pack_weight_kernel<<<blocks, 256, 0,
                         at::cuda::getCurrentCUDAStream()>>>(
        weight.data_ptr<float>(),
        reinterpret_cast<half*>(packed.data_ptr<at::Half>()),
        vectors);
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return packed;
}

struct MatmulPlan {
    cublasLtMatmulDesc_t operation = nullptr;
    cublasLtMatrixLayout_t a = nullptr;
    cublasLtMatrixLayout_t b = nullptr;
    cublasLtMatrixLayout_t c = nullptr;
    cublasLtMatmulAlgo_t algorithm{};
    size_t workspace_bytes = 0;

    ~MatmulPlan() {
        if (operation) cublasLtMatmulDescDestroy(operation);
        if (a) cublasLtMatrixLayoutDestroy(a);
        if (b) cublasLtMatrixLayoutDestroy(b);
        if (c) cublasLtMatrixLayoutDestroy(c);
    }
};

struct ThreadContext {
    cublasLtHandle_t handle = nullptr;
    std::map<std::tuple<int, int, int, int>,
             std::unique_ptr<MatmulPlan>> plans;

    ThreadContext() {
        check_lt(cublasLtCreate(&handle));
    }

    ~ThreadContext() {
        plans.clear();
        if (handle) cublasLtDestroy(handle);
    }

    MatmulPlan& get_plan(
        int device, int m, int n, int k, const float* bias
    ) {
        auto key = std::make_tuple(device, m, n, k);
        auto found = plans.find(key);
        if (found != plans.end())
            return *found->second;

        auto plan = std::make_unique<MatmulPlan>();
        check_lt(cublasLtMatmulDescCreate(
            &plan->operation, CUBLAS_COMPUTE_32F, CUDA_R_32F));

        cublasOperation_t trans_a = CUBLAS_OP_T;
        cublasOperation_t trans_b = CUBLAS_OP_N;
        cublasLtEpilogue_t epilogue = CUBLASLT_EPILOGUE_RELU_BIAS;

        check_lt(cublasLtMatmulDescSetAttribute(
            plan->operation, CUBLASLT_MATMUL_DESC_TRANSA,
            &trans_a, sizeof(trans_a)));
        check_lt(cublasLtMatmulDescSetAttribute(
            plan->operation, CUBLASLT_MATMUL_DESC_TRANSB,
            &trans_b, sizeof(trans_b)));
        check_lt(cublasLtMatmulDescSetAttribute(
            plan->operation, CUBLASLT_MATMUL_DESC_EPILOGUE,
            &epilogue, sizeof(epilogue)));
        check_lt(cublasLtMatmulDescSetAttribute(
            plan->operation, CUBLASLT_MATMUL_DESC_BIAS_POINTER,
            &bias, sizeof(bias)));

        // Column-major interpretation gives output = weight @ input.T.
        check_lt(cublasLtMatrixLayoutCreate(
            &plan->a, CUDA_R_16F, k, n, k));
        check_lt(cublasLtMatrixLayoutCreate(
            &plan->b, CUDA_R_16F, k, m, k));
        check_lt(cublasLtMatrixLayoutCreate(
            &plan->c, CUDA_R_32F, n, m, n));

        cublasLtMatmulPreference_t preference;
        check_lt(cublasLtMatmulPreferenceCreate(&preference));
        const size_t workspace_limit = 64ULL << 20;
        auto status = cublasLtMatmulPreferenceSetAttribute(
            preference, CUBLASLT_MATMUL_PREF_MAX_WORKSPACE_BYTES,
            &workspace_limit, sizeof(workspace_limit));
        if (status != CUBLAS_STATUS_SUCCESS) {
            cublasLtMatmulPreferenceDestroy(preference);
            check_lt(status);
        }

        cublasLtMatmulHeuristicResult_t result{};
        int count = 0;
        status = cublasLtMatmulAlgoGetHeuristic(
            handle, plan->operation,
            plan->a, plan->b, plan->c, plan->c,
            preference, 1, &result, &count);
        cublasLtMatmulPreferenceDestroy(preference);
        check_lt(status);
        TORCH_CHECK(count > 0, "No compatible cuBLASLt algorithm");
        check_lt(result.state);

        plan->algorithm = result.algo;
        plan->workspace_bytes = result.workspaceSize;

        auto inserted = plans.emplace(key, std::move(plan));
        return *inserted.first->second;
    }
};

torch::Tensor linear_activate_cuda(
    torch::Tensor x,
    torch::Tensor packed_weight,
    torch::Tensor bias,
    double subtract_value,
    double multiply_value
) {
    const c10::cuda::CUDAGuard guard(x.device());

    // clone also handles contiguous tensors with unaligned storage offsets.
    if (!x.is_contiguous() ||
        (reinterpret_cast<uintptr_t>(x.data_ptr()) & 15))
        x = x.clone(at::MemoryFormat::Contiguous);
    if (!bias.is_contiguous() ||
        (reinterpret_cast<uintptr_t>(bias.data_ptr()) & 15))
        bias = bias.clone(at::MemoryFormat::Contiguous);

    const int m = static_cast<int>(x.size(0));
    const int n = static_cast<int>(packed_weight.size(0));
    const int k = static_cast<int>(x.size(1));

    auto output = torch::empty({m, n}, x.options());
    auto xh = torch::empty(x.sizes(), x.options().dtype(at::kHalf));
    auto adjusted_bias = torch::empty({n}, bias.options());

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    const int64_t vectors = x.numel() / 4;
    const int64_t work = std::max<int64_t>(vectors, n / 4);
    const int blocks = static_cast<int>(
        std::min<int64_t>((work + 255) / 256, 4096));

    pack_input_bias_kernel<<<blocks, 256, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<half*>(xh.data_ptr<at::Half>()),
        bias.data_ptr<float>(),
        adjusted_bias.data_ptr<float>(),
        vectors, n / 4,
        static_cast<float>(subtract_value),
        static_cast<float>(multiply_value));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    static thread_local ThreadContext context;
    const float* bias_pointer = adjusted_bias.data_ptr<float>();
    MatmulPlan& plan = context.get_plan(
        x.get_device(), m, n, k, bias_pointer);
    check_lt(cublasLtMatmulDescSetAttribute(
        plan.operation, CUBLASLT_MATMUL_DESC_BIAS_POINTER,
        &bias_pointer, sizeof(bias_pointer)));

    auto workspace = torch::empty(
        {static_cast<int64_t>(plan.workspace_bytes)},
        x.options().dtype(at::kByte));

    const float alpha = static_cast<float>(multiply_value) / 256.0f;
    const float beta = 0.0f;
    check_lt(cublasLtMatmul(
        context.handle, plan.operation,
        &alpha,
        packed_weight.data_ptr<at::Half>(), plan.a,
        xh.data_ptr<at::Half>(), plan.b,
        &beta,
        output.data_ptr<float>(), plan.c,
        output.data_ptr<float>(), plan.c,
        &plan.algorithm,
        workspace.data_ptr(), plan.workspace_bytes,
        stream));

    return output;
}
"""

_extension = load_inline(
    name="linear_packed_lt_fused_v5",
    cpp_sources=r"""
torch::Tensor pack_weight_cuda(torch::Tensor weight);
torch::Tensor linear_activate_cuda(
    torch::Tensor x,
    torch::Tensor packed_weight,
    torch::Tensor bias,
    double subtract_value,
    double multiply_value);
""",
    cuda_sources=_cuda_source,
    functions=["pack_weight_cuda", "linear_activate_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    extra_ldflags=["-lcublasLt"],
    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
        self.register_buffer("_packed_weight", None, persistent=False)
        self._packed_weight_key = None

    def _apply(self, fn, recurse=True):
        self._packed_weight = None
        self._packed_weight_key = None
        return super()._apply(fn, recurse=recurse)

    def forward(self, x):
        weight = self.linear.weight
        bias = self.linear.bias
        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or weight.dtype != torch.float32
            or bias.dtype != torch.float32
            or torch.is_autocast_enabled()
            or self.linear.in_features == 0
            or self.linear.out_features == 0
            or self.linear.in_features % 8 != 0
            or self.linear.out_features % 8 != 0
            or x.numel() == 0
        ):
            return torch.relu(
                (F.linear(x, weight, bias) - self.subtract_value)
                * self.multiply_value
            )

        if x.ndim == 0 or x.shape[-1] != self.linear.in_features:
            raise ValueError("Input's last dimension must match in_features")
        if x.device != weight.device or x.device != bias.device:
            raise ValueError("Input and parameters must share a device")

        try:
            version = weight._version
        except RuntimeError:
            version = None

        stream = torch.cuda.current_stream(weight.device)
        key = (
            id(weight),
            weight.data_ptr(),
            version,
            tuple(weight.shape),
            tuple(weight.stride()),
            weight.device,
            stream.cuda_stream,
        )

        if (
            version is None
            or self._packed_weight is None
            or self._packed_weight_key != key
        ):
            source = weight
            if not source.is_contiguous() or source.data_ptr() % 16:
                source = source.clone(memory_format=torch.contiguous_format)
            packed_weight = _extension.pack_weight_cuda(source)
            if version is not None:
                self._packed_weight = packed_weight
                self._packed_weight_key = key
        else:
            packed_weight = self._packed_weight

        shape = (*x.shape[:-1], self.linear.out_features)
        rows = math.prod(x.shape[:-1])
        return _extension.linear_activate_cuda(
            x.reshape(rows, self.linear.in_features),
            packed_weight,
            bias,
            float(self.subtract_value),
            float(self.multiply_value),
        ).view(shape)
```