```python
# Copyright (c) Meta Platforms, Inc. and affiliates.
# All rights reserved.
#
# This source code is licensed under the license found in the
# LICENSE file in the root directory of this source tree.

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

# Fused CUDA kernel:
#   out = hardtanh(gelu(tanh(swish(h + add_value))))
# The GEMM (x @ W^T + bias) is executed with cuBLAS via at::addmm, and the
# entire element-wise chain (bias-add of add_value, Swish, Tanh, GELU,
# Hardtanh) is fused into a single kernel, eliminating ~5 extra global-memory
# round trips and several kernel launches.
cuda_source = r"""
#include <torch/extension.h>
#include <ATen/ATen.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda_runtime.h>

__device__ __forceinline__ float fused_act(float v) {
    // Swish: v * sigmoid(v) = v / (1 + exp(-v))
    float sv = v / (1.0f + expf(-v));
    // Tanh
    float tv = tanhf(sv);
    // GELU (exact erf formulation, matches F.gelu default)
    float gv = 0.5f * tv * (1.0f + erff(tv * 0.70710678118654752440f));
    // Hardtanh clamp to [-1, 1]
    return fminf(fmaxf(gv, -1.0f), 1.0f);
}

__global__ void fused_epilogue_scalar_kernel(
    const float* __restrict__ h,
    const float* __restrict__ addv,
    float* __restrict__ out,
    int cols)
{
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y;
    if (col < cols) {
        long idx = (long)row * cols + col;
        out[idx] = fused_act(h[idx] + addv[col]);
    }
}

__global__ void fused_epilogue_vec4_kernel(
    const float4* __restrict__ h,
    const float4* __restrict__ addv,
    float4* __restrict__ out,
    int cols4)
{
    int c4 = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y;
    if (c4 < cols4) {
        long idx = (long)row * cols4 + c4;
        float4 val = h[idx];
        float4 av  = addv[c4];
        val.x = fused_act(val.x + av.x);
        val.y = fused_act(val.y + av.y);
        val.z = fused_act(val.z + av.z);
        val.w = fused_act(val.w + av.w);
        out[idx] = val;
    }
}

torch::Tensor fused_forward_cuda(torch::Tensor x, torch::Tensor weight, torch::Tensor bias, torch::Tensor add_value) {
    TORCH_CHECK(x.is_cuda() && weight.is_cuda() && bias.is_cuda() && add_value.is_cuda(),
                "all inputs must be CUDA tensors");
    TORCH_CHECK(x.scalar_type() == at::kFloat, "only float32 is supported");

    // GEMM with bias via cuBLAS (same numerics as nn.Linear).
    auto h = at::addmm(bias, x, weight.t());  // (batch, out_features)
    auto out = torch::empty_like(h);

    int rows = h.size(0);
    int cols = h.size(1);
    if (rows == 0 || cols == 0) {
        return out;
    }

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

    if ((cols % 4) == 0) {
        int cols4 = cols / 4;
        dim3 grid((cols4 + threads - 1) / threads, rows);
        fused_epilogue_vec4_kernel<<<grid, threads, 0, stream>>>(
            reinterpret_cast<const float4*>(h.data_ptr<float>()),
            reinterpret_cast<const float4*>(add_value.data_ptr<float>()),
            reinterpret_cast<float4*>(out.data_ptr<float>()),
            cols4);
    } else {
        dim3 grid((cols + threads - 1) / threads, rows);
        fused_epilogue_scalar_kernel<<<grid, threads, 0, stream>>>(
            h.data_ptr<float>(),
            add_value.data_ptr<float>(),
            out.data_ptr<float>(),
            cols);
    }
    return out;
}
"""

cpp_source = (
    "torch::Tensor fused_forward_cuda(torch::Tensor x, torch::Tensor weight, "
    "torch::Tensor bias, torch::Tensor add_value);"
)

fused_module = load_inline(
    name="fused_linear_swish_tanh_gelu_hardtanh",
    cpp_sources=cpp_source,
    cuda_sources=cuda_source,
    functions=["fused_forward_cuda"],
    verbose=False,
    extra_cuda_cflags=["-O3"],
)


class ModelNew(nn.Module):
    def __init__(self, in_features, out_features, add_value_shape):
        super(ModelNew, self).__init__()
        self.matmul = nn.Linear(in_features, out_features)
        self.add_value = nn.Parameter(torch.randn(add_value_shape))
        self.fused = fused_module

    def forward(self, x):
        if x.is_cuda and x.dtype == torch.float32:
            return self.fused.fused_forward_cuda(
                x, self.matmul.weight, self.matmul.bias, self.add_value
            )
        # Fallback (e.g. CPU tensors): original operator sequence.
        x = F.linear(x, self.matmul.weight, self.matmul.bias)
        x = x + self.add_value
        x = torch.sigmoid(x) * x  # Swish
        x = torch.tanh(x)
        x = F.gelu(x)  # GELU
        x = F.hardtanh(x, min_val=-1, max_val=1)  # Hardtanh
        return x
```