I want to quickly solve linear least squares problem for $x \in \mathbb{R}^n$

$$min_x \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 michg 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.