```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 pool_sum_scale_cuda(
    torch::Tensor input,
    int64_t kernel_size,
    double scale_factor);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/AccumulateType.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <limits>

template <typename T>
__device__ __forceinline__ T warp_sum(T value) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        value += __shfl_down_sync(0xffffffff, value, offset);
    }
    return value;
}

template <typename scalar_t, int FixedKernel>
__global__ void pool_sum_scale_kernel(
    const scalar_t* __restrict__ input,
    scalar_t* __restrict__ output,
    int64_t features,
    int64_t pooled_features,
    int64_t runtime_kernel,
    at::acc_type<scalar_t, true> scale)
{
    using acc_t = at::acc_type<scalar_t, true>;
    constexpr int Threads = 256;
    __shared__ acc_t warp_sums[Threads / 32];

    const int tid = threadIdx.x;
    const int64_t kernel =
        FixedKernel > 0 ? FixedKernel : runtime_kernel;
    const scalar_t* row = input + int64_t(blockIdx.x) * features;

    acc_t total = acc_t(0);

    for (int64_t p = tid; p < pooled_features; p += Threads) {
        const int64_t start = p * kernel;
        acc_t maximum = -std::numeric_limits<acc_t>::infinity();

        if (FixedKernel == 2) {
            const acc_t a = static_cast<acc_t>(row[start]);
            const acc_t b = static_cast<acc_t>(row[start + 1]);
            maximum = (a > b || isnan(a)) ? a : b;
        } else {
            for (int64_t j = 0; j < kernel; ++j) {
                const acc_t value = static_cast<acc_t>(row[start + j]);
                if (value > maximum || isnan(value)) {
                    maximum = value;
                }
            }
        }
        total += maximum;
    }

    total = warp_sum(total);

    if ((tid & 31) == 0) {
        warp_sums[tid >> 5] = total;
    }
    __syncthreads();

    if (tid < 32) {
        acc_t result = tid < Threads / 32 ? warp_sums[tid] : acc_t(0);
        result = warp_sum(result);

        if (tid == 0) {
            // Preserve the intermediate dtype of torch.sum before scaling.
            const scalar_t rounded_sum = static_cast<scalar_t>(result);
            output[blockIdx.x] = static_cast<scalar_t>(
                static_cast<acc_t>(rounded_sum) * scale);
        }
    }
}

torch::Tensor pool_sum_scale_cuda(
    torch::Tensor input,
    int64_t kernel_size,
    double scale_factor)
{
    TORCH_CHECK(input.is_cuda(), "input must be a CUDA tensor");
    TORCH_CHECK(input.dim() == 2, "input must have shape (batch, features)");
    TORCH_CHECK(input.is_contiguous(), "input must be contiguous");
    TORCH_CHECK(kernel_size > 0, "kernel_size must be positive");
    TORCH_CHECK(input.size(1) >= kernel_size, "pooling output is empty");

    const c10::cuda::CUDAGuard device_guard(input.device());
    const int64_t batch = input.size(0);
    const int64_t features = input.size(1);
    const int64_t pooled_features = features / kernel_size;

    auto output = torch::empty({batch}, input.options());
    if (batch == 0) {
        return output;
    }

    const auto stream = at::cuda::getCurrentCUDAStream();

    AT_DISPATCH_FLOATING_TYPES_AND2(
        at::ScalarType::Half,
        at::ScalarType::BFloat16,
        input.scalar_type(),
        "pool_sum_scale_cuda",
        [&] {
            using acc_t = at::acc_type<scalar_t, true>;
            const acc_t scale = static_cast<acc_t>(scale_factor);

            if (kernel_size == 2) {
                pool_sum_scale_kernel<scalar_t, 2>
                    <<<batch, 256, 0, stream>>>(
                        input.data_ptr<scalar_t>(),
                        output.data_ptr<scalar_t>(),
                        features,
                        pooled_features,
                        kernel_size,
                        scale);
            } else {
                pool_sum_scale_kernel<scalar_t, 0>
                    <<<batch, 256, 0, stream>>>(
                        input.data_ptr<scalar_t>(),
                        output.data_ptr<scalar_t>(),
                        features,
                        pooled_features,
                        kernel_size,
                        scale);
            }
        });

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="fused_maxpool_sum_scale_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["pool_sum_scale_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self, in_features, out_features, kernel_size, scale_factor):
        super().__init__()
        self.matmul = nn.Linear(in_features, out_features)
        self.kernel_size = int(kernel_size)
        self.scale_factor = float(scale_factor)

    def forward(self, x):
        x = self.matmul(x)
        return _fused_ops.pool_sum_scale_cuda(
            x, self.kernel_size, self.scale_factor
        )
```