```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 <array>
#include <cstdint>
#include <limits>
#include <map>
#include <memory>
#include <tuple>
#include <vector>

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 count
) {
    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < count;
         i += int64_t(blockDim.x) * gridDim.x) {
        const float4 v = reinterpret_cast<const float4*>(src)[i];
        half2* p = reinterpret_cast<half2*>(dst) + 2 * i;
        p[0] = __floats2half2_rn(v.x * 256.0f, v.y * 256.0f);
        p[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 count,
    int bias_count,
    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 < count; i += step) {
        const float4 v = reinterpret_cast<const float4*>(src)[i];
        half2* p = reinterpret_cast<half2*>(dst) + 2 * i;
        p[0] = __floats2half2_rn(v.x, v.y);
        p[1] = __floats2half2_rn(v.z, v.w);
    }

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

__global__ void finish_kernel(
    float* __restrict__ output,
    const float* __restrict__ bias,
    int64_t count,
    int columns
) {
    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < count;
         i += int64_t(blockDim.x) * gridDim.x) {
        float4 v = reinterpret_cast<float4*>(output)[i];
        const float4 b =
            reinterpret_cast<const float4*>(bias)[i % columns];
        v.x = __fadd_rn(v.x, b.x);
        v.y = __fadd_rn(v.y, b.y);
        v.z = __fadd_rn(v.z, b.z);
        v.w = __fadd_rn(v.w, b.w);
        v.x = v.x < 0.0f ? 0.0f : v.x;
        v.y = v.y < 0.0f ? 0.0f : v.y;
        v.z = v.z < 0.0f ? 0.0f : v.z;
        v.w = v.w < 0.0f ? 0.0f : v.w;
        reinterpret_cast<float4*>(output)[i] = v;
    }
}

torch::Tensor pack_weight_cuda(torch::Tensor weight) {
    const c10::cuda::CUDAGuard guard(weight.device());
    auto packed = torch::empty(
        weight.sizes(), weight.options().dtype(at::kHalf));

    const int64_t count = weight.numel() / 4;
    const int blocks = static_cast<int>(
        std::min<int64_t>((count + 255) / 256, 8192));

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

struct Candidate {
    cublasLtMatmulAlgo_t algorithm{};
    size_t workspace_bytes = 0;
    int operation = 0;
};

struct MatmulPlan {
    std::array<cublasLtMatmulDesc_t, 2> operations{{nullptr, nullptr}};
    cublasLtMatrixLayout_t a = nullptr;
    cublasLtMatrixLayout_t b = nullptr;
    cublasLtMatrixLayout_t c = nullptr;
    std::vector<Candidate> candidates;
    size_t workspace_bytes = 0;
    int selected = 0;
    bool tuned = false;

    ~MatmulPlan() {
        for (auto operation : operations)
            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
    ) {
        const 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(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;
        check_lt(cublasLtMatmulPreferenceSetAttribute(
            preference, CUBLASLT_MATMUL_PREF_MAX_WORKSPACE_BYTES,
            &workspace_limit, sizeof(workspace_limit)));

        for (int mode = 0; mode < 2; ++mode) {
            auto& operation = plan->operations[mode];
            check_lt(cublasLtMatmulDescCreate(
                &operation, CUBLAS_COMPUTE_32F, CUDA_R_32F));

            cublasOperation_t trans_a = CUBLAS_OP_T;
            cublasOperation_t trans_b = CUBLAS_OP_N;
            check_lt(cublasLtMatmulDescSetAttribute(
                operation, CUBLASLT_MATMUL_DESC_TRANSA,
                &trans_a, sizeof(trans_a)));
            check_lt(cublasLtMatmulDescSetAttribute(
                operation, CUBLASLT_MATMUL_DESC_TRANSB,
                &trans_b, sizeof(trans_b)));

            if (mode == 0) {
                cublasLtEpilogue_t epilogue =
                    CUBLASLT_EPILOGUE_RELU_BIAS;
                check_lt(cublasLtMatmulDescSetAttribute(
                    operation, CUBLASLT_MATMUL_DESC_EPILOGUE,
                    &epilogue, sizeof(epilogue)));
                check_lt(cublasLtMatmulDescSetAttribute(
                    operation, CUBLASLT_MATMUL_DESC_BIAS_POINTER,
                    &bias, sizeof(bias)));
            }

            std::array<cublasLtMatmulHeuristicResult_t, 16> results{};
            int count = 0;
            const auto status = cublasLtMatmulAlgoGetHeuristic(
                handle, operation,
                plan->a, plan->b, plan->c, plan->c,
                preference, static_cast<int>(results.size()),
                results.data(), &count);

            if (status != CUBLAS_STATUS_SUCCESS)
                continue;

            for (int i = 0; i < count; ++i) {
                if (results[i].state != CUBLAS_STATUS_SUCCESS)
                    continue;
                Candidate candidate;
                candidate.algorithm = results[i].algo;
                candidate.workspace_bytes = results[i].workspaceSize;
                candidate.operation = mode;
                plan->candidates.push_back(candidate);
                plan->workspace_bytes = std::max(
                    plan->workspace_bytes, candidate.workspace_bytes);
            }
        }

        cublasLtMatmulPreferenceDestroy(preference);
        TORCH_CHECK(!plan->candidates.empty(),
                    "No compatible cuBLASLt algorithm");

        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());

    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 input_vectors = x.numel() / 4;
    const int64_t work = std::max<int64_t>(input_vectors, n / 4);
    const int pack_blocks = static_cast<int>(
        std::min<int64_t>((work + 255) / 256, 4096));

    pack_input_bias_kernel<<<pack_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>(),
        input_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.operations[0], CUBLASLT_MATMUL_DESC_BIAS_POINTER,
        &bias_pointer, sizeof(bias_pointer)));

    const size_t workspace_bytes = plan.tuned
        ? plan.candidates[plan.selected].workspace_bytes
        : plan.workspace_bytes;
    auto workspace = torch::empty(
        {static_cast<int64_t>(workspace_bytes)},
        x.options().dtype(at::kByte));

    const float alpha = static_cast<float>(multiply_value) / 256.0f;
    const float beta = 0.0f;
    const int64_t output_vectors = int64_t(m) * n / 4;
    const int finish_blocks = static_cast<int>(
        std::min<int64_t>((output_vectors + 255) / 256, 4096));

    auto run = [&](const Candidate& candidate) {
        const auto status = cublasLtMatmul(
            context.handle, plan.operations[candidate.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,
            &candidate.algorithm,
            workspace.data_ptr(), candidate.workspace_bytes,
            stream);

        if (status == CUBLAS_STATUS_SUCCESS &&
            candidate.operation == 1) {
            finish_kernel<<<finish_blocks, 256, 0, stream>>>(
                output.data_ptr<float>(), bias_pointer,
                output_vectors, n / 4);
        }
        return status;
    };

    cudaStreamCaptureStatus capture_status;
    C10_CUDA_CHECK(cudaStreamIsCapturing(stream, &capture_status));

    if (!plan.tuned && capture_status == cudaStreamCaptureStatusNone) {
        cudaEvent_t start, stop;
        C10_CUDA_CHECK(cudaEventCreate(&start));
        C10_CUDA_CHECK(cudaEventCreate(&stop));

        float best_time = std::numeric_limits<float>::infinity();
        int best = -1;

        for (int i = 0; i < static_cast<int>(plan.candidates.size()); ++i) {
            const auto& candidate = plan.candidates[i];
            if (run(candidate) != CUBLAS_STATUS_SUCCESS)
                continue;

            C10_CUDA_CHECK(cudaEventRecord(start, stream));
            bool valid = true;
            for (int repeat = 0; repeat < 3; ++repeat) {
                if (run(candidate) != CUBLAS_STATUS_SUCCESS) {
                    valid = false;
                    break;
                }
            }
            C10_CUDA_CHECK(cudaEventRecord(stop, stream));
            C10_CUDA_CHECK(cudaEventSynchronize(stop));

            float elapsed = 0.0f;
            C10_CUDA_CHECK(cudaEventElapsedTime(&elapsed, start, stop));
            if (valid && elapsed < best_time) {
                best_time = elapsed;
                best = i;
            }
        }

        C10_CUDA_CHECK(cudaEventDestroy(start));
        C10_CUDA_CHECK(cudaEventDestroy(stop));
        TORCH_CHECK(best >= 0, "No executable cuBLASLt algorithm");
        plan.selected = best;
        plan.tuned = true;
    }

    check_lt(run(plan.candidates[plan.selected]));
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_extension = load_inline(
    name="linear_packed_lt_adaptive_v6",
    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

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