We introduce an algorithm for reliably computing quantities associated with several types of semiparametric mixed models in situations where the condition number on the random effects matrix is large. The algorithm is numerically stable and efficient. It was designed to process penalized spline (P-spline) models without making unnecessary numerical approximations. The algorithm, PSQR (P-splines via QR), is formulated in terms of QR decompositions. PSQR can treat both exactly rank deficient and ill-conditioned matrices. The latter situation often arises in large scale mixed models and/or when a P-spline is estimated using a basis with poor numerical properties, e.g. a truncated power function (TPF) basis. We provide concrete examples where unnecessary numerical approximations introduce both subtle and dramatic errors that would likely go undetected, thus demonstrating the importance of using this reliable numerical algorithm. Simulation results studying a univariate function and a longitudinal data set are used to demonstrate the algorithm. Extensions and the utility of the method in more general semiparametric regression applications are briefly discussed. MATLAB scripts demonstrating implementation are provided in the Supplemental Materials.