We present a novel variant of the least squares method for approximating elliptic partial differential equations of second order using hierarchical tensor product splines on two- and three-dimensional grids. It is based on collocation of both the PDE and the boundary conditions using appropriate weighting. Hence, neither meshing nor numerical integration is necessary. Using C2-splines, we can compute the residual pointwise. On one hand, these values are convenient for steering adaptive refinement. On the other hand, they admit the a posteriori estimation of the approximation error via the maximum principle. Special care is taken to make sure that the so-derived error bounds are not unduly pessimistic. As a result, we obtain an efficient algorithm that is able to compute approximate solutions satisfying prescribed bounds on the error with respect to the maximum norm.