```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_forward_cuda(
    torch::Tensor input,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    int64_t groups,
    double eps);
"""

_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 <cublas_v2.h>

template<bool INPUT>
__global__ void split_pack_kernel(
    const float* __restrict__ source,
    __half* __restrict__ packed,
    int64_t count,
    int features)
{
    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < count;
         i += int64_t(blockDim.x) * gridDim.x) {
        const float value = source[i];
        const __half high = __float2half_rn(value);
        const __half low =
            __float2half_rn((value - __half2float(high)) * 4096.0f);

        const int64_t row = i / features;
        const int col = i - row * features;
        const int64_t offset = row * (2LL * features) + col;

        if (INPUT) {
            packed[offset] = low;
            packed[offset + features] = high;
        } else {
            packed[offset] = high;
            packed[offset + features] = low;
        }
    }
}

__device__ __forceinline__ float warp_sum(float value) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        value += __shfl_down_sync(0xffffffff, value, offset);
    return __shfl_sync(0xffffffff, value, 0);
}

__device__ __forceinline__ float finish(
    float value,
    float mean,
    float inv_std,
    float gamma,
    float beta,
    float multiplier)
{
    float normalized = ((value - mean) * inv_std) * gamma + beta;
    float activated = normalized * (1.0f / (1.0f + expf(-normalized)));
    float weighted = activated * multiplier;
    return weighted * (1.0f / (1.0f + expf(-weighted)));
}

template<bool SIZE_32>
__global__ void norm_activation_kernel(
    const float* __restrict__ input,
    const float* __restrict__ bias,
    const float* __restrict__ gamma,
    const float* __restrict__ beta,
    const float* __restrict__ multiplier,
    float* __restrict__ output,
    int channels,
    int groups,
    int group_size,
    int total_groups,
    float eps)
{
    const int lane = threadIdx.x & 31;
    const int group = blockIdx.x * (blockDim.x / 32) + threadIdx.x / 32;
    if (group >= total_groups) return;

    const int channel_start = (group % groups) * group_size;
    const int64_t base =
        int64_t(group / groups) * channels + channel_start;

    if (SIZE_32) {
        const int channel = channel_start + lane;
        const float value = input[base + lane] + bias[channel];
        const float mean = warp_sum(value) * 0.03125f;
        const float centered = value - mean;
        const float variance = warp_sum(centered * centered) * 0.03125f;
        const float inv_std = rsqrtf(variance + eps);

        output[base + lane] = finish(
            value, mean, inv_std,
            gamma[channel], beta[channel], multiplier[channel]);
    } else {
        float sum = 0.0f;
        for (int i = lane; i < group_size; i += 32)
            sum += input[base + i] + bias[channel_start + i];

        const float mean = warp_sum(sum) / float(group_size);
        float squared_sum = 0.0f;
        for (int i = lane; i < group_size; i += 32) {
            const float centered =
                input[base + i] + bias[channel_start + i] - mean;
            squared_sum += centered * centered;
        }
        const float inv_std =
            rsqrtf(warp_sum(squared_sum) / float(group_size) + eps);

        for (int i = lane; i < group_size; i += 32) {
            const int channel = channel_start + i;
            output[base + i] = finish(
                input[base + i] + bias[channel],
                mean, inv_std,
                gamma[channel], beta[channel], multiplier[channel]);
        }
    }
}

torch::Tensor fused_forward_cuda(
    torch::Tensor input,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    int64_t groups,
    double eps)
{
    const c10::cuda::CUDAGuard guard(input.device());

    TORCH_CHECK(input.is_cuda() && input.scalar_type() == at::kFloat);
    TORCH_CHECK(input.dim() == 2 && input.is_contiguous());
    TORCH_CHECK(weight.dim() == 2);

    const int batch = static_cast<int>(input.size(0));
    const int features = static_cast<int>(input.size(1));
    const int channels = static_cast<int>(weight.size(0));

    TORCH_CHECK(features > 0 && weight.size(1) == features);
    TORCH_CHECK(groups > 0 && channels % groups == 0);

    for (const auto& tensor : {weight, bias, gamma, beta, multiplier}) {
        TORCH_CHECK(tensor.device() == input.device());
        TORCH_CHECK(tensor.scalar_type() == at::kFloat);
        TORCH_CHECK(tensor.is_contiguous());
    }
    TORCH_CHECK(bias.numel() == channels);
    TORCH_CHECK(gamma.numel() == channels);
    TORCH_CHECK(beta.numel() == channels);
    TORCH_CHECK(multiplier.numel() == channels);

    auto output = torch::empty({batch, channels}, input.options());
    if (batch == 0) return output;

    auto intermediate = torch::empty_like(output);
    auto half_options = input.options().dtype(at::kHalf);
    auto packed_input = torch::empty({batch, 2LL * features}, half_options);
    auto packed_weight = torch::empty({channels, 2LL * features}, half_options);

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

    const int64_t input_count = input.numel();
    const int64_t weight_count = weight.numel();
    const int input_blocks =
        static_cast<int>(std::min<int64_t>((input_count + 255) / 256, 8192));
    const int weight_blocks =
        static_cast<int>(std::min<int64_t>((weight_count + 255) / 256, 8192));

    auto* xp = reinterpret_cast<__half*>(packed_input.data_ptr<at::Half>());
    auto* wp = reinterpret_cast<__half*>(packed_weight.data_ptr<at::Half>());

    split_pack_kernel<true><<<input_blocks, threads, 0, stream>>>(
        input.data_ptr<float>(), xp, input_count, features);
    split_pack_kernel<false><<<weight_blocks, threads, 0, stream>>>(
        weight.data_ptr<float>(), wp, weight_count, features);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    auto handle = at::cuda::getCurrentCUDABlasHandle();

    cublasMath_t original_math;
    TORCH_CHECK(cublasGetMathMode(handle, &original_math) ==
                CUBLAS_STATUS_SUCCESS);

    const auto precise_math = static_cast<cublasMath_t>(
        static_cast<int>(original_math) |
        static_cast<int>(CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION));

    TORCH_CHECK(cublasSetMathMode(handle, precise_math) ==
                CUBLAS_STATUS_SUCCESS);

    const float correction_scale = 1.0f / 4096.0f;
    const float zero = 0.0f;
    const float one = 1.0f;
    const int leading_dim = 2 * features;

    // Packed GEMM computes both cross terms:
    // W_high @ X_low.T + W_low @ X_high.T.
    const cublasStatus_t correction_status = cublasGemmEx(
        handle,
        CUBLAS_OP_T, CUBLAS_OP_N,
        channels, batch, 2 * features,
        &correction_scale,
        wp, CUDA_R_16F, leading_dim,
        xp, CUDA_R_16F, leading_dim,
        &zero,
        intermediate.data_ptr<float>(), CUDA_R_32F, channels,
        CUBLAS_COMPUTE_32F,
        CUBLAS_GEMM_DEFAULT_TENSOR_OP);

    // Accumulate the dominant product in FP32.
    cublasStatus_t main_status = correction_status;
    if (correction_status == CUBLAS_STATUS_SUCCESS) {
        main_status = cublasGemmEx(
            handle,
            CUBLAS_OP_T, CUBLAS_OP_N,
            channels, batch, features,
            &one,
            wp, CUDA_R_16F, leading_dim,
            xp + features, CUDA_R_16F, leading_dim,
            &one,
            intermediate.data_ptr<float>(), CUDA_R_32F, channels,
            CUBLAS_COMPUTE_32F,
            CUBLAS_GEMM_DEFAULT_TENSOR_OP);
    }

    const cublasStatus_t restore_status =
        cublasSetMathMode(handle, original_math);

    TORCH_CHECK(correction_status == CUBLAS_STATUS_SUCCESS,
                "Correction GEMM failed: ", int(correction_status));
    TORCH_CHECK(main_status == CUBLAS_STATUS_SUCCESS,
                "Main GEMM failed: ", int(main_status));
    TORCH_CHECK(restore_status == CUBLAS_STATUS_SUCCESS,
                "Could not restore cuBLAS math mode");

    const int group_size = channels / static_cast<int>(groups);
    const int total_groups = batch * static_cast<int>(groups);
    const int blocks = (total_groups + 7) / 8;

    if (group_size == 32) {
        norm_activation_kernel<true><<<blocks, threads, 0, stream>>>(
            intermediate.data_ptr<float>(),
            bias.data_ptr<float>(),
            gamma.data_ptr<float>(),
            beta.data_ptr<float>(),
            multiplier.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, static_cast<int>(groups), group_size,
            total_groups, static_cast<float>(eps));
    } else {
        norm_activation_kernel<false><<<blocks, threads, 0, stream>>>(
            intermediate.data_ptr<float>(),
            bias.data_ptr<float>(),
            gamma.data_ptr<float>(),
            beta.data_ptr<float>(),
            multiplier.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, static_cast<int>(groups), group_size,
            total_groups, static_cast<float>(eps));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="compensated_half_gemm_fused_gn_swish_v3",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_forward_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    extra_ldflags=["-lcublas"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_features, out_features, num_groups, multiply_weight_shape
    ):
        super().__init__()
        self.gemm = nn.Linear(in_features, out_features)
        self.group_norm = nn.GroupNorm(num_groups, out_features)
        self.multiply_weight = nn.Parameter(torch.randn(multiply_weight_shape))

    def forward(self, x):
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 2
            and x.shape[1] > 0
            and x.shape[1] % 8 == 0
            and self.gemm.out_features % 8 == 0
            and self.gemm.weight.dtype == torch.float32
            and self.gemm.bias.dtype == torch.float32
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
            and self.multiply_weight.dtype == torch.float32
            and self.multiply_weight.shape == (self.gemm.out_features,)
            and not torch.is_autocast_enabled()
        ):
            return _fused_ops.fused_forward_cuda(
                x.contiguous(),
                self.gemm.weight.contiguous(),
                self.gemm.bias.contiguous(),
                self.group_norm.weight.contiguous(),
                self.group_norm.bias.contiguous(),
                self.multiply_weight.contiguous(),
                self.group_norm.num_groups,
                self.group_norm.eps,
            )

        x = self.gemm(x)
        x = self.group_norm(x)
        x = x * torch.sigmoid(x)
        x = x * self.multiply_weight
        return x * torch.sigmoid(x)
```