This is just a large scale optimization problem with a quadratic objective function and very simple inequality constraints. Your problem can be written as
$\min$ $x^T(A^TA)x -x^TA^Tb-b^TAx + b^Tb$ subject to $u-x\geq0$ and $x-l\geq0$
This is such a well-solved problem that there are many software packages that have built in functions for this. Generally it is not worth thinking about writing your own algorithm because experts have already spent vasted amount of time on this problem and have written algorithms which will out-perform anything you can write in a few days.
Here are a selection of built-in functions in different software packages that can be used:
In Matlab: quadprog lsqlin (type help quadprog lsqlin into Matlab and it will tell tells you exactly what options are best to use.type. I have just (approximately) solved your problem with random sparse matrices and it works great.)
KNITRO for Mathematica this package also solves this exact problem but I don't have this software so I can't tell you which exact function.
For a free solver I have found this: http://sourceforge.net/projects/quadprog/ However it assumes that $A$ has full column rankso . This is just because this algorithm uses the dual problem which exists when the Hessian $A^TA$ is positive definite(this ensures the dual problem exists) then a free script written in C++ can be downloaded from here: http://sourceforge.net/projects/quadprog/.