```python
import torch
import torch.nn as nn
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(0xffffffff, value, offset);
    }
    return value;
}

template<bool PAIRS>
__global__ void bias_pool_sum_scale_kernel(
    const float* __restrict__ input,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int features,
    int kernel_size,
    float scale
) {
    const int row = blockIdx.x;
    const int tid = threadIdx.x;
    const int pooled_features = features / kernel_size;
    const float* row_input =
        input + static_cast<int64_t>(row) * features;

    float sum = 0.0f;

    if (PAIRS) {
        const float2* values =
            reinterpret_cast<const float2*>(row_input);
        const float2* biases =
            reinterpret_cast<const float2*>(bias);

        for (int p = tid; p < pooled_features; p += blockDim.x) {
            const float2 v = values[p];
            const float2 b = biases[p];
            const float first = v.x + b.x;
            const float second = v.y + b.y;
            sum += (isnan(first) || first > second) ? first : second;
        }
    } else {
        for (int p = tid; p < pooled_features; p += blockDim.x) {
            float maximum = -CUDART_INF_F;
            const int start = p * kernel_size;
            for (int j = 0; j < kernel_size; ++j) {
                const float value = row_input[start + j] + bias[start + j];
                if (isnan(value) || value > maximum) {
                    maximum = value;
                }
            }
            sum += maximum;
        }
    }

    sum = warp_sum(sum);

    __shared__ float warp_sums[8];
    const int lane = tid & 31;
    const int warp = tid >> 5;

    if (lane == 0) {
        warp_sums[warp] = sum;
    }
    __syncthreads();

    if (warp == 0) {
        float total = lane < 8 ? warp_sums[lane] : 0.0f;
        total = warp_sum(total);
        if (lane == 0) {
            output[row] = total * scale;
        }
    }
}

torch::Tensor bias_pool_sum_scale_cuda(
    torch::Tensor input,
    torch::Tensor bias,
    int64_t kernel_size,
    double scale
) {
    TORCH_CHECK(input.is_cuda() && bias.is_cuda(),
                "Expected CUDA tensors");
    TORCH_CHECK(input.device() == bias.device(),
                "Input and bias must be on the same device");
    TORCH_CHECK(input.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat,
                "Expected float32 tensors");
    TORCH_CHECK(input.dim() == 2 && bias.dim() == 1,
                "Expected a matrix and a bias vector");
    TORCH_CHECK(input.is_contiguous() && bias.is_contiguous(),
                "Expected contiguous tensors");
    TORCH_CHECK(input.size(1) == bias.numel(),
                "Bias size mismatch");
    TORCH_CHECK(kernel_size > 0 && kernel_size <= input.size(1),
                "Invalid pooling kernel size");
    TORCH_CHECK(input.size(1) <= 2147483647,
                "Feature dimension is too large");

    const c10::cuda::CUDAGuard guard(input.device());
    auto output = torch::empty({input.size(0)}, input.options());

    if (input.size(0) == 0) {
        return output;
    }

    const int features = static_cast<int>(input.size(1));
    const int pool_size = static_cast<int>(kernel_size);
    const dim3 grid(static_cast<unsigned int>(input.size(0)));
    const auto stream = c10::cuda::getCurrentCUDAStream().stream();

    if (pool_size == 2 && features % 2 == 0) {
        bias_pool_sum_scale_kernel<true><<<grid, 256, 0, stream>>>(
            input.data_ptr<float>(),
            bias.data_ptr<float>(),
            output.data_ptr<float>(),
            features,
            pool_size,
            static_cast<float>(scale)
        );
    } else {
        bias_pool_sum_scale_kernel<false><<<grid, 256, 0, stream>>>(
            input.data_ptr<float>(),
            bias.data_ptr<float>(),
            output.data_ptr<float>(),
            features,
            pool_size,
            static_cast<float>(scale)
        );
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor bias_pool_sum_scale_cuda(
    torch::Tensor input,
    torch::Tensor bias,
    int64_t kernel_size,
    double scale
);
"""

_fused_ops = load_inline(
    name="linear_bias_pool_sum_scale_ext",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["bias_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):
        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or self.matmul.weight.dtype != torch.float32
            or torch.is_autocast_enabled()
        ):
            y = self.matmul(x)
            y = self.max_pool(y.unsqueeze(1)).squeeze(1)
            return y.sum(dim=1) * self.scale_factor

        # Retain cuBLAS GEMM; fuse bias, pooling, reduction, and scaling.
        y = torch.mm(x, self.matmul.weight.t())
        return _fused_ops.bias_pool_sum_scale_cuda(
            y,
            self.matmul.bias,
            self.kernel_size,
            self.scale_factor,
        )
```