```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 x,
    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>

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

template <typename scalar_t>
__device__ __forceinline__ scalar_t finish(
    float value,
    scalar_t gamma,
    scalar_t beta,
    scalar_t multiplier,
    float mean,
    float inv_std)
{
    // Preserve the intermediate dtype of the original PyTorch operations.
    scalar_t normalized = static_cast<scalar_t>(
        (value - mean) * inv_std * static_cast<float>(gamma)
        + static_cast<float>(beta));

    scalar_t sigmoid1 = static_cast<scalar_t>(
        1.0f / (1.0f + __expf(-static_cast<float>(normalized))));

    scalar_t swish1 = static_cast<scalar_t>(
        static_cast<float>(normalized) * static_cast<float>(sigmoid1));

    scalar_t product = static_cast<scalar_t>(
        static_cast<float>(swish1) * static_cast<float>(multiplier));

    scalar_t sigmoid2 = static_cast<scalar_t>(
        1.0f / (1.0f + __expf(-static_cast<float>(product))));

    return static_cast<scalar_t>(
        static_cast<float>(product) * static_cast<float>(sigmoid2));
}

template <typename scalar_t, bool GROUP_SIZE_32>
__global__ void fused_kernel(
    const scalar_t* __restrict__ x,
    const scalar_t* __restrict__ gamma,
    const scalar_t* __restrict__ beta,
    const scalar_t* __restrict__ multiplier,
    scalar_t* __restrict__ output,
    int64_t total_groups,
    int groups_per_row,
    int group_size,
    float eps)
{
    const int lane = threadIdx.x & 31;
    const int64_t group =
        static_cast<int64_t>(blockIdx.x) * (blockDim.x / 32)
        + threadIdx.x / 32;

    if (group >= total_groups)
        return;

    const int64_t base = group * group_size;
    const int channel_base = (group % groups_per_row) * group_size;

    if (GROUP_SIZE_32) {
        const float value = static_cast<float>(x[base + lane]);
        const float mean = warp_sum(value) * (1.0f / 32.0f);
        const float centered = value - mean;
        const float variance =
            warp_sum(centered * centered) * (1.0f / 32.0f);
        const float inv_std = rsqrtf(variance + eps);
        const int channel = channel_base + lane;

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

        const float mean = warp_sum(sum) / group_size;

        float squared_sum = 0.0f;
        for (int i = lane; i < group_size; i += 32) {
            const float centered = static_cast<float>(x[base + i]) - mean;
            squared_sum += centered * centered;
        }

        const float inv_std =
            rsqrtf(warp_sum(squared_sum) / group_size + eps);

        for (int i = lane; i < group_size; i += 32) {
            const int channel = channel_base + i;
            output[base + i] = finish(
                static_cast<float>(x[base + i]),
                gamma[channel], beta[channel], multiplier[channel],
                mean, inv_std);
        }
    }
}

torch::Tensor fused_groupnorm_swish_cuda(
    torch::Tensor x,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    int64_t groups,
    double eps)
{
    TORCH_CHECK(x.is_cuda() && x.is_contiguous() && x.dim() == 2,
                "x must be a contiguous two-dimensional CUDA tensor");
    TORCH_CHECK(groups > 0 && x.size(1) % groups == 0,
                "channels must be divisible by groups");

    const int64_t channels = x.size(1);
    for (const auto& parameter : {gamma, beta, multiplier}) {
        TORCH_CHECK(parameter.device() == x.device(),
                    "all tensors must be on the same device");
        TORCH_CHECK(parameter.scalar_type() == x.scalar_type(),
                    "all tensors must have the same dtype");
        TORCH_CHECK(parameter.is_contiguous() &&
                    parameter.numel() == channels,
                    "parameters must be contiguous channel vectors");
    }

    const c10::cuda::CUDAGuard device_guard(x.device());
    auto output = torch::empty_like(x);
    const int64_t total_groups = x.size(0) * groups;
    if (total_groups == 0)
        return output;

    constexpr int threads = 256;
    const int blocks = static_cast<int>((total_groups + 7) / 8);
    const int group_size = static_cast<int>(channels / groups);
    const auto stream = at::cuda::getCurrentCUDAStream();

    AT_DISPATCH_FLOATING_TYPES_AND2(
        at::ScalarType::Half,
        at::ScalarType::BFloat16,
        x.scalar_type(),
        "fused_groupnorm_swish",
        [&] {
            if (group_size == 32) {
                fused_kernel<scalar_t, true><<<blocks, threads, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    gamma.data_ptr<scalar_t>(),
                    beta.data_ptr<scalar_t>(),
                    multiplier.data_ptr<scalar_t>(),
                    output.data_ptr<scalar_t>(),
                    total_groups, static_cast<int>(groups), group_size,
                    static_cast<float>(eps));
            } else {
                fused_kernel<scalar_t, false><<<blocks, threads, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    gamma.data_ptr<scalar_t>(),
                    beta.data_ptr<scalar_t>(),
                    multiplier.data_ptr<scalar_t>(),
                    output.data_ptr<scalar_t>(),
                    total_groups, static_cast<int>(groups), group_size,
                    static_cast<float>(eps));
            }
        });

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="fused_groupnorm_double_swish",
    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):
        x = self.gemm(x)

        parameters = (
            self.group_norm.weight,
            self.group_norm.bias,
            self.multiply_weight,
        )
        if (
            x.is_cuda
            and x.dim() == 2
            and x.is_contiguous()
            and x.dtype in (torch.float32, torch.float16, torch.bfloat16)
            and self.multiply_weight.shape == (x.shape[1],)
            and all(
                p.dtype == x.dtype
                and p.device == x.device
                and p.is_contiguous()
                for p in parameters
            )
        ):
            return _fused_ops.fused_groupnorm_swish_cuda(
                x,
                self.group_norm.weight,
                self.group_norm.bias,
                self.multiply_weight,
                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)
```