I want to quickly solve the following linear least-squares problem

$$\min_{x \in \mathbb{R}^n} \left\| A x - b \right\|_2^2$$

with a special sparse structure where each row in $A$ has only up to $4$ consecutive non-zero entries. This makes its normal matrix

$$C = A^T A$$

a positive definite $7$-*band* matrix with a condition number between $8^2$ and $400^2$. So, these condition numbers aren't too bad that solving the Gauß normal equation system

$$ C x = A^T b $$

instead wouldn't get me into trouble numerically, I think. What are my options? I could try conjugate gradient methods but I would prefer direct solvers that can deal with these kinds of special cases in $O(n)$ time independent of the condition. I'm aware of algorithms for the tridiagonal case and I guess I could try to adapt them for $5$ bands (?). But before reinventing the wheel and/or testing many different algorithms I wanted to ask you about what approach might be the most efficient in terms of time because I have lots of these problems (millions) with values for $n$ of around $5000$ where $A$ has about $4n$ rows.