```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);
"""

_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.h>
#include <cuda_runtime.h>
#include <math.h>

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, typename acc_t, int StaticK>
__global__ void pool_sum_scale_kernel(
    const scalar_t* __restrict__ input,
    scalar_t* __restrict__ output,
    int64_t width,
    int64_t pooled_width,
    int64_t dynamic_k,
    acc_t scale)
{
    constexpr int Threads = 256;
    __shared__ acc_t warp_totals[Threads / 32];

    const int tid = threadIdx.x;
    const int64_t k = StaticK ? StaticK : dynamic_k;
    const scalar_t* row = input + int64_t(blockIdx.x) * width;

    acc_t sums[4] = {acc_t(0), acc_t(0), acc_t(0), acc_t(0)};

    for (int64_t base = tid; base < pooled_width; base += Threads * 4) {
        #pragma unroll
        for (int lane = 0; lane < 4; ++lane) {
            const int64_t p = base + lane * Threads;
            if (p < pooled_width) {
                acc_t maximum = -acc_t(INFINITY);

                for (int64_t j = 0; j < k; ++j) {
                    const acc_t value = static_cast<acc_t>(row[p * k + j]);
                    if (isnan(value) || value > maximum)
                        maximum = value;
                }

                sums[lane] += maximum;
            }
        }
    }

    acc_t total = (sums[0] + sums[1]) + (sums[2] + sums[3]);
    total = warp_sum(total);

    if ((tid & 31) == 0)
        warp_totals[tid >> 5] = total;
    __syncthreads();

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

        if (tid == 0) {
            // torch.sum returns the input dtype before the scaling operation.
            const scalar_t reduced = static_cast<scalar_t>(total);
            output[blockIdx.x] =
                static_cast<scalar_t>(static_cast<acc_t>(reduced) * scale);
        }
    }
}

torch::Tensor pool_sum_scale_cuda(
    torch::Tensor input, int64_t kernel_size, double scale)
{
    TORCH_CHECK(input.is_cuda(), "input must be a CUDA tensor");
    TORCH_CHECK(input.dim() == 2, "input must have two dimensions");
    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 rows = input.size(0);
    const int64_t width = input.size(1);
    const int64_t pooled_width = width / kernel_size;

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

    const cudaStream_t stream =
        at::cuda::getCurrentCUDAStream(input.get_device());

    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>;

            if (kernel_size == 2) {
                pool_sum_scale_kernel<scalar_t, acc_t, 2>
                    <<<rows, 256, 0, stream>>>(
                        input.data_ptr<scalar_t>(),
                        output.data_ptr<scalar_t>(),
                        width,
                        pooled_width,
                        kernel_size,
                        static_cast<acc_t>(scale));
            } else {
                pool_sum_scale_kernel<scalar_t, acc_t, 0>
                    <<<rows, 256, 0, stream>>>(
                        input.data_ptr<scalar_t>(),
                        output.data_ptr<scalar_t>(),
                        width,
                        pooled_width,
                        kernel_size,
                        static_cast<acc_t>(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.max_pool = nn.MaxPool1d(kernel_size)
        self.kernel_size = int(kernel_size)
        self.scale_factor = float(scale_factor)

    def forward(self, x):
        x = self.matmul(x)
        if not x.is_cuda:
            x = self.max_pool(x.unsqueeze(1)).squeeze(1)
            return x.sum(dim=1) * self.scale_factor

        return _fused_ops.pool_sum_scale_cuda(
            x, self.kernel_size, self.scale_factor
        )
```