```python
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 <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAStream.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);
}

__device__ __forceinline__ float swish(float value) {
    return value / (1.0f + __expf(-value));
}

template<int StaticGroupSize>
__global__ void fused_groupnorm_kernel(
    const float* __restrict__ input,
    const float* __restrict__ gamma,
    const float* __restrict__ beta,
    const float* __restrict__ multiplier,
    float* __restrict__ output,
    int64_t total_groups,
    int num_groups,
    int dynamic_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 int group_size =
        StaticGroupSize > 0 ? StaticGroupSize : dynamic_group_size;
    const int64_t base = group * group_size;
    const int channel_base = (group % num_groups) * group_size;

    if (StaticGroupSize == 32) {
        const float value = input[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 normalized = centered * rsqrtf(variance + eps);
        const int channel = channel_base + lane;
        const float affine = normalized * gamma[channel] + beta[channel];
        const float scaled = swish(affine) * multiplier[channel];
        output[base + lane] = swish(scaled);
    } else {
        float sum = 0.0f;
        for (int i = lane; i < group_size; i += 32) {
            sum += input[base + i];
        }
        const float mean = warp_sum(sum) / static_cast<float>(group_size);

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

        for (int i = lane; i < group_size; i += 32) {
            const int channel = channel_base + i;
            const float normalized = (input[base + i] - mean) * inv_std;
            const float affine =
                normalized * gamma[channel] + beta[channel];
            const float scaled = swish(affine) * multiplier[channel];
            output[base + i] = swish(scaled);
        }
    }
}

torch::Tensor fused_groupnorm_swish(
    torch::Tensor input,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    int64_t num_groups,
    double eps
) {
    TORCH_CHECK(input.is_cuda(), "input must be a CUDA tensor");
    TORCH_CHECK(input.dim() == 2, "input must have shape (N, C)");
    TORCH_CHECK(input.scalar_type() == at::kFloat,
                "input must have float32 dtype");
    TORCH_CHECK(input.is_contiguous(), "input must be contiguous");

    const int64_t channels = input.size(1);
    TORCH_CHECK(num_groups > 0 && channels % num_groups == 0,
                "channels must be divisible by num_groups");

    for (const auto& parameter : {gamma, beta, multiplier}) {
        TORCH_CHECK(parameter.device() == input.device(),
                    "parameters must be on the input device");
        TORCH_CHECK(parameter.scalar_type() == at::kFloat,
                    "parameters must have float32 dtype");
        TORCH_CHECK(parameter.is_contiguous(),
                    "parameters must be contiguous");
        TORCH_CHECK(parameter.numel() == channels,
                    "parameters must contain C elements");
    }

    const c10::cuda::CUDAGuard device_guard(input.device());
    auto output = torch::empty_like(input);
    const int64_t total_groups = input.size(0) * num_groups;
    if (total_groups == 0) {
        return output;
    }

    constexpr int threads = 256;
    constexpr int warps_per_block = threads / 32;
    const int blocks =
        static_cast<int>((total_groups + warps_per_block - 1)
                         / warps_per_block);
    const int group_size = static_cast<int>(channels / num_groups);
    const auto stream = c10::cuda::getCurrentCUDAStream();

    if (group_size == 32) {
        fused_groupnorm_kernel<32><<<blocks, threads, 0, stream.stream()>>>(
            input.data_ptr<float>(),
            gamma.data_ptr<float>(),
            beta.data_ptr<float>(),
            multiplier.data_ptr<float>(),
            output.data_ptr<float>(),
            total_groups,
            static_cast<int>(num_groups),
            group_size,
            static_cast<float>(eps)
        );
    } else {
        fused_groupnorm_kernel<0><<<blocks, threads, 0, stream.stream()>>>(
            input.data_ptr<float>(),
            gamma.data_ptr<float>(),
            beta.data_ptr<float>(),
            multiplier.data_ptr<float>(),
            output.data_ptr<float>(),
            total_groups,
            static_cast<int>(num_groups),
            group_size,
            static_cast<float>(eps)
        );
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_cpp_source = r"""
torch::Tensor fused_groupnorm_swish(
    torch::Tensor input,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    int64_t num_groups,
    double eps
);
"""

_fused_ops = load_inline(
    name="gemm_groupnorm_double_swish_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_groupnorm_swish"],
    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)

        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 == (x.shape[1],)
        ):
            return _fused_ops.fused_groupnorm_swish(
                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)
```