```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 multiply_weight,
    int64_t groups,
    double eps);
"""

_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 activate(
    float value, float gamma, float beta, float multiplier,
    float mean, float inv_std) {
    float z = (value - mean) * inv_std;
    z = z * gamma + beta;
    z = z / (1.0f + expf(-z));
    z *= multiplier;
    return z / (1.0f + expf(-z));
}

template<bool GROUP_SIZE_32>
__global__ void fused_groupnorm_swish_kernel(
    const float* __restrict__ x,
    const float* __restrict__ gamma,
    const float* __restrict__ beta,
    const float* __restrict__ multiplier,
    float* __restrict__ output,
    int64_t total_groups,
    int channels,
    int groups,
    int group_size,
    int64_t multiplier_stride0,
    int64_t multiplier_stride1,
    float eps) {

    const int lane = threadIdx.x & 31;
    const int64_t group_index =
        static_cast<int64_t>(blockIdx.x) * (blockDim.x / 32)
        + threadIdx.x / 32;

    // Each warp owns one complete normalization group.
    if (group_index >= total_groups)
        return;

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

    if (GROUP_SIZE_32) {
        const int channel = channel_start + lane;
        const float value = 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 float m = multiplier[
            row * multiplier_stride0 + channel * multiplier_stride1];

        output[base + lane] = activate(
            value, gamma[channel], beta[channel], m, mean, inv_std);
    } else {
        float sum = 0.0f;
        for (int j = lane; j < group_size; j += 32)
            sum += x[base + j];

        const float mean = warp_sum(sum) / group_size;

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

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

        for (int j = lane; j < group_size; j += 32) {
            const int channel = channel_start + j;
            const float m = multiplier[
                row * multiplier_stride0 + channel * multiplier_stride1];

            output[base + j] = activate(
                x[base + j], gamma[channel], beta[channel],
                m, mean, inv_std);
        }
    }
}

torch::Tensor fused_groupnorm_swish_cuda(
    torch::Tensor x,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiply_weight,
    int64_t groups,
    double eps) {

    TORCH_CHECK(x.is_cuda() && x.scalar_type() == torch::kFloat32,
                "x must be a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 2 && x.is_contiguous(),
                "x must be contiguous and two-dimensional");

    const c10::cuda::CUDAGuard device_guard(x.device());
    const int channels = static_cast<int>(x.size(1));

    TORCH_CHECK(groups > 0 && channels % groups == 0,
                "channels must be divisible by groups");
    TORCH_CHECK(
        gamma.device() == x.device() &&
        beta.device() == x.device() &&
        multiply_weight.device() == x.device(),
        "all tensors must be on the same device");
    TORCH_CHECK(
        gamma.scalar_type() == torch::kFloat32 &&
        beta.scalar_type() == torch::kFloat32 &&
        multiply_weight.scalar_type() == torch::kFloat32,
        "parameters must have float32 dtype");
    TORCH_CHECK(
        gamma.is_contiguous() && beta.is_contiguous() &&
        gamma.numel() == channels && beta.numel() == channels,
        "invalid GroupNorm parameters");

    auto multiplier = multiply_weight.expand_as(x);
    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;
    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 = channels / groups;
    const auto stream = c10::cuda::getCurrentCUDAStream();

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

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="fused_groupnorm_double_swish_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_groupnorm_swish_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3", "--fmad=false"],
    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 self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
            and self.multiply_weight.dtype == torch.float32
        ):
            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)
```