```python
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.cpp_extension import load_inline

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

std::vector<torch::Tensor> linear_pool_sum_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t kernel_size,
    double scale);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <math_constants.h>
#include <climits>
#include <vector>

static void check_cublas(cublasStatus_t status) {
    TORCH_CHECK(status == CUBLAS_STATUS_SUCCESS,
                "cuBLAS failed with status ", static_cast<int>(status));
}

__device__ __forceinline__ float maximum_nan(float a, float b) {
    return isnan(b) || a < b ? b : a;
}

__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_kernel(
    const float* __restrict__ input,
    const float* __restrict__ bias,
    float* __restrict__ output,
    int columns,
    int kernel_size,
    float scale
) {
    const int tid = threadIdx.x;
    const float* row = input + static_cast<int64_t>(blockIdx.x) * columns;

    float sum0 = 0.0f;
    float sum1 = 0.0f;

    if constexpr (PAIRS) {
        const float4* values = reinterpret_cast<const float4*>(row);
        const float4* biases = reinterpret_cast<const float4*>(bias);
        for (int i = tid; i < columns / 4; i += blockDim.x) {
            const float4 v = values[i];
            const float4 b = biases[i];
            sum0 += maximum_nan(v.x + b.x, v.y + b.y);
            sum1 += maximum_nan(v.z + b.z, v.w + b.w);
        }
    } else {
        for (int i = tid; i < columns / kernel_size; i += blockDim.x) {
            const int start = i * kernel_size;
            float maximum = -CUDART_INF_F;
            for (int j = 0; j < kernel_size; ++j)
                maximum = maximum_nan(
                    maximum, row[start + j] + bias[start + j]);
            sum0 += maximum;
        }
    }

    float sum = warp_sum(sum0 + sum1);
    __shared__ float partial[8];

    const int lane = tid & 31;
    const int warp = tid >> 5;
    if (lane == 0)
        partial[warp] = sum;
    __syncthreads();

    if (warp == 0) {
        sum = warp_sum(lane < 8 ? partial[lane] : 0.0f);
        if (lane == 0)
            output[blockIdx.x] = sum * scale;
    }
}

std::vector<torch::Tensor> linear_pool_sum_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t kernel_size,
    double scale
) {
    TORCH_CHECK(x.is_cuda() && weight.is_cuda() && bias.is_cuda(),
                "All inputs must be CUDA tensors");
    TORCH_CHECK(x.device() == weight.device() &&
                x.device() == bias.device(), "Device mismatch");
    TORCH_CHECK(x.scalar_type() == at::kFloat &&
                weight.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat,
                "All inputs must be float32");
    TORCH_CHECK(x.dim() == 2 && weight.dim() == 2 && bias.dim() == 1,
                "Invalid input dimensions");
    TORCH_CHECK(x.is_contiguous() && weight.is_contiguous() &&
                bias.is_contiguous(), "Inputs must be contiguous");
    TORCH_CHECK(x.size(1) == weight.size(1) &&
                bias.size(0) == weight.size(0), "Shape mismatch");
    TORCH_CHECK(kernel_size > 0 && kernel_size <= weight.size(0),
                "Invalid pooling kernel size");
    TORCH_CHECK(x.size(0) <= INT_MAX && x.size(1) <= INT_MAX &&
                weight.size(0) <= INT_MAX, "Dimensions are too large");

    const c10::cuda::CUDAGuard guard(x.device());
    const int batch = static_cast<int>(x.size(0));
    const int inputs = static_cast<int>(x.size(1));
    const int outputs = static_cast<int>(weight.size(0));

    auto linear = torch::empty({batch, outputs}, x.options());
    auto result = torch::empty({batch}, x.options());

    if (batch == 0)
        return {result, linear};

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

    if (inputs == 0) {
        C10_CUDA_CHECK(cudaMemsetAsync(
            linear.data_ptr<float>(), 0,
            linear.numel() * sizeof(float), stream));
    } else {
        cublasHandle_t handle = at::cuda::getCurrentCUDABlasHandle();
        cublasMath_t previous_math;
        cublasPointerMode_t previous_pointer;
        check_cublas(cublasGetMathMode(handle, &previous_math));
        check_cublas(cublasGetPointerMode(handle, &previous_pointer));
        check_cublas(cublasSetPointerMode(handle, CUBLAS_POINTER_MODE_HOST));
        check_cublas(cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH));

        const float alpha = 1.0f;
        const float beta = 0.0f;

        // Row-major X @ W.T is column-major W @ X.T.
        const cublasStatus_t status = cublasGemmEx(
            handle,
            CUBLAS_OP_T, CUBLAS_OP_N,
            outputs, batch, inputs,
            &alpha,
            weight.data_ptr<float>(), CUDA_R_32F, inputs,
            x.data_ptr<float>(), CUDA_R_32F, inputs,
            &beta,
            linear.data_ptr<float>(), CUDA_R_32F, outputs,
            CUBLAS_COMPUTE_32F_FAST_TF32,
            CUBLAS_GEMM_DEFAULT_TENSOR_OP);

        const auto math_status = cublasSetMathMode(handle, previous_math);
        const auto pointer_status =
            cublasSetPointerMode(handle, previous_pointer);
        check_cublas(status);
        check_cublas(math_status);
        check_cublas(pointer_status);
    }

    if (kernel_size == 2 && outputs % 4 == 0 &&
        reinterpret_cast<uintptr_t>(bias.data_ptr<float>()) % 16 == 0) {
        bias_pool_sum_kernel<true><<<batch, 256, 0, stream>>>(
            linear.data_ptr<float>(), bias.data_ptr<float>(),
            result.data_ptr<float>(), outputs, 2,
            static_cast<float>(scale));
    } else {
        bias_pool_sum_kernel<false><<<batch, 256, 0, stream>>>(
            linear.data_ptr<float>(), bias.data_ptr<float>(),
            result.data_ptr<float>(), outputs,
            static_cast<int>(kernel_size), static_cast<float>(scale));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return {result, linear};
}
"""

_ops = load_inline(
    name="tf32_linear_fused_bias_pool_sum_v3",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["linear_pool_sum_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    extra_ldflags=["-lcublas"],
    verbose=False,
)


class _LinearPoolSum(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, kernel_size, scale):
        result, linear = _ops.linear_pool_sum_cuda(
            x, weight, bias, kernel_size, scale
        )
        ctx.kernel_size = kernel_size
        ctx.scale = scale
        if any(ctx.needs_input_grad[:3]):
            ctx.save_for_backward(x, weight, bias, linear)
        return result

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, bias, linear = ctx.saved_tensors
        _, indices = F.max_pool1d(
            (linear + bias).unsqueeze(1),
            kernel_size=ctx.kernel_size,
            stride=ctx.kernel_size,
            return_indices=True,
        )
        indices = indices.squeeze(1)
        grad_linear = torch.zeros_like(linear)
        values = (grad_output * ctx.scale).unsqueeze(1).expand_as(indices)
        grad_linear.scatter_(1, indices, values)

        grad_x = (
            grad_linear.matmul(weight) if ctx.needs_input_grad[0] else None
        )
        grad_weight = (
            grad_linear.t().matmul(x) if ctx.needs_input_grad[1] else None
        )
        grad_bias = (
            grad_linear.sum(dim=0) if ctx.needs_input_grad[2] else None
        )
        return grad_x, grad_weight, grad_bias, None, None


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 = (
            kernel_size[0]
            if isinstance(kernel_size, (tuple, list))
            else kernel_size
        )
        self.scale_factor = scale_factor

    def forward(self, x):
        weight = self.matmul.weight
        bias = self.matmul.bias
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 2
            and weight.dtype == torch.float32
            and bias.dtype == torch.float32
            and not torch.is_autocast_enabled()
        ):
            return _LinearPoolSum.apply(
                x.contiguous(),
                weight.contiguous(),
                bias.contiguous(),
                self.kernel_size,
                self.scale_factor,
            )

        x = self.matmul(x)
        x = self.max_pool(x.unsqueeze(1)).squeeze(1)
        return x.sum(dim=1) * self.scale_factor
```