```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,
    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 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 activate(
    float centered,
    float inv_std,
    float gamma,
    float beta,
    float multiplier)
{
    const float normalized = (centered * inv_std) * gamma + beta;
    const float swish = normalized / (1.0f + __expf(-normalized));
    const float weighted = swish * multiplier;
    return weighted / (1.0f + __expf(-weighted));
}

__global__ void fused_groupnorm32_kernel(
    const float4* __restrict__ input,
    const float4* __restrict__ gamma,
    const float4* __restrict__ beta,
    const float4* __restrict__ multiplier,
    float4* __restrict__ output,
    int64_t vector_count,
    int channel_vectors,
    float eps)
{
    const int64_t index =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;

    // Each normalization group occupies eight consecutive threads.
    // Tensor dimensions guarantee that any partial warp contains whole groups.
    const bool valid = index < vector_count;
    float4 value = make_float4(0.f, 0.f, 0.f, 0.f);
    if (valid)
        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 local_variance =
        (value.x * value.x + value.y * value.y) +
        (value.z * value.z + value.w * value.w);

    const float inv_std =
        rsqrtf(sum8(local_variance) * 0.03125f + eps);

    if (valid) {
        const int channel = static_cast<int>(index % channel_vectors);
        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;
    }
}

torch::Tensor fused_groupnorm_swish_cuda(
    torch::Tensor input,
    torch::Tensor gamma,
    torch::Tensor beta,
    torch::Tensor multiplier,
    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());
    const int64_t channels = input.size(1);
    TORCH_CHECK(channels > 0 && channels % 32 == 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());
        TORCH_CHECK(tensor.numel() == channels);
    }

    auto output = torch::empty_like(input);
    const int64_t vector_count = input.numel() / 4;
    if (vector_count == 0)
        return output;

    constexpr int threads = 256;
    const int blocks =
        static_cast<int>((vector_count + threads - 1) / threads);
    const auto stream = at::cuda::getCurrentCUDAStream();

    fused_groupnorm32_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>()),
        vector_count,
        static_cast<int>(channels / 4),
        static_cast<float>(eps));

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="groupnorm32_double_swish_vectorized_v4",
    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):
        # Retain the optimized GEMM and its configured precision policy.
        x = self.gemm(x)

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 2
            and x.is_contiguous()
            and self.gemm.out_features // self.group_norm.num_groups == 32
            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,
                self.group_norm.weight.contiguous(),
                self.group_norm.bias.contiguous(),
                self.multiply_weight.contiguous(),
                self.group_norm.eps,
            )

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