```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_groupnorm_swish_cuda(
    torch::Tensor input,
    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 <algorithm>

__device__ __forceinline__ float sum8(float value) {
    value += __shfl_xor_sync(0xffffffffu, value, 4, 8);
    value += __shfl_xor_sync(0xffffffffu, value, 2, 8);
    value += __shfl_xor_sync(0xffffffffu, value, 1, 8);
    return value;
}

__device__ __forceinline__ float sum32(float value) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        value += __shfl_xor_sync(0xffffffffu, value, offset);
    return value;
}

__device__ __forceinline__ float activate(
    float centered,
    float inv_std,
    float gamma,
    float beta,
    float multiplier)
{
    const float normalized = (centered * inv_std) * gamma + beta;
    const float first = normalized * (1.0f / (1.0f + __expf(-normalized)));
    const float weighted = first * multiplier;
    return weighted * (1.0f / (1.0f + __expf(-weighted)));
}

__global__ void groupnorm32_swish_kernel(
    const float4* __restrict__ input,
    const float4* __restrict__ gamma,
    const float4* __restrict__ beta,
    const float4* __restrict__ multiplier,
    float4* __restrict__ output,
    int64_t vector_count,
    int64_t channel_vectors,
    float eps)
{
    for (int64_t index = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         index < vector_count;
         index += int64_t(blockDim.x) * gridDim.x) {
        // Each eight-thread subgroup owns one complete group of 32 values.
        const int64_t channel = index % channel_vectors;
        float4 value = input[index];

        const float mean =
            sum8((value.x + value.y) + (value.z + value.w)) * 0.03125f;

        value.x -= mean;
        value.y -= mean;
        value.z -= mean;
        value.w -= mean;

        const float squares =
            (value.x * value.x + value.y * value.y) +
            (value.z * value.z + value.w * value.w);
        const float inv_std = rsqrtf(sum8(squares) * 0.03125f + eps);

        const float4 g = gamma[channel];
        const float4 b = beta[channel];
        const float4 m = multiplier[channel];

        float4 result;
        result.x = activate(value.x, inv_std, g.x, b.x, m.x);
        result.y = activate(value.y, inv_std, g.y, b.y, m.y);
        result.z = activate(value.z, inv_std, g.z, b.z, m.z);
        result.w = activate(value.w, inv_std, g.w, b.w, m.w);
        output[index] = result;
    }
}

__global__ void groupnorm_swish_kernel(
    const float* __restrict__ input,
    const float* __restrict__ gamma,
    const float* __restrict__ beta,
    const float* __restrict__ multiplier,
    float* __restrict__ output,
    int64_t groups,
    int64_t group_size,
    int64_t total_groups,
    float eps)
{
    const int lane = threadIdx.x & 31;
    const int64_t first_group =
        int64_t(blockIdx.x) * (blockDim.x / 32) + threadIdx.x / 32;
    const int64_t group_stride = int64_t(gridDim.x) * (blockDim.x / 32);

    for (int64_t group = first_group;
         group < total_groups;
         group += group_stride) {
        const int64_t base = group * group_size;
        const int64_t channel_base = (group % groups) * group_size;

        float sum = 0.0f;
        for (int64_t i = lane; i < group_size; i += 32)
            sum += input[base + i];

        const float mean = sum32(sum) / float(group_size);

        float squares = 0.0f;
        for (int64_t i = lane; i < group_size; i += 32) {
            const float centered = input[base + i] - mean;
            squares += centered * centered;
        }

        const float inv_std =
            rsqrtf(sum32(squares) / float(group_size) + eps);

        for (int64_t i = lane; i < group_size; i += 32) {
            const int64_t channel = channel_base + i;
            output[base + i] = activate(
                input[base + i] - mean,
                inv_std,
                gamma[channel],
                beta[channel],
                multiplier[channel]);
        }
    }
}

torch::Tensor fused_groupnorm_swish_cuda(
    torch::Tensor input,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    int64_t groups,
    double eps)
{
    TORCH_CHECK(input.is_cuda() && input.scalar_type() == at::kFloat);
    TORCH_CHECK(input.dim() == 2 && input.is_contiguous());
    const c10::cuda::CUDAGuard guard(input.device());

    const int64_t channels = input.size(1);
    TORCH_CHECK(groups > 0 && channels > 0 && channels % groups == 0);

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

    auto output = torch::empty_like(input);
    if (input.numel() == 0)
        return output;

    const auto stream = at::cuda::getCurrentCUDAStream();
    constexpr int threads = 256;
    const int64_t group_size = channels / groups;

    const bool aligned =
        (reinterpret_cast<uintptr_t>(input.data_ptr<float>()) % 16 == 0) &&
        (reinterpret_cast<uintptr_t>(gamma.data_ptr<float>()) % 16 == 0) &&
        (reinterpret_cast<uintptr_t>(beta.data_ptr<float>()) % 16 == 0) &&
        (reinterpret_cast<uintptr_t>(multiplier.data_ptr<float>()) % 16 == 0);

    if (group_size == 32 && aligned) {
        const int64_t count = input.numel() / 4;
        const int blocks =
            int(std::min<int64_t>((count + threads - 1) / threads, 8192));

        groupnorm32_swish_kernel<<<blocks, threads, 0, stream>>>(
            reinterpret_cast<const float4*>(input.data_ptr<float>()),
            reinterpret_cast<const float4*>(gamma.data_ptr<float>()),
            reinterpret_cast<const float4*>(beta.data_ptr<float>()),
            reinterpret_cast<const float4*>(multiplier.data_ptr<float>()),
            reinterpret_cast<float4*>(output.data_ptr<float>()),
            count,
            channels / 4,
            static_cast<float>(eps));
    } else {
        const int64_t total_groups = input.size(0) * groups;
        const int blocks =
            int(std::min<int64_t>((total_groups + 7) / 8, 8192));

        groupnorm_swish_kernel<<<blocks, threads, 0, stream>>>(
            input.data_ptr<float>(),
            gamma.data_ptr<float>(),
            beta.data_ptr<float>(),
            multiplier.data_ptr<float>(),
            output.data_ptr<float>(),
            groups,
            group_size,
            total_groups,
            static_cast<float>(eps));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="fused_groupnorm_double_swish_v6",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_groupnorm_swish_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    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):
        # Preserve the reference GEMM's precision and optimized backend.
        x = self.gemm(x)

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 2
            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,)
        ):
            return _fused_ops.fused_groupnorm_swish_cuda(
                x.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.group_norm(x)
        x = x * torch.sigmoid(x)
        x = x * self.multiply_weight
        return x * torch.sigmoid(x)
```