The computation of a few singular triplets of large, sparse matrices is a challenging task, especially when the smallest magnitude singular values are needed in high accuracy. Most recent efforts try to address this problem through variations of the Lanczos bidiagonalization method, but algorithmic research is ongoing and without production level software. We show that a more efficient, robust, and production-ready method can be developed on top of the state-of-the-art eigensolver PRIMME. In a first stage, our new method, primme_svds, solves the normal equations problem up to the best achievable accuracy. If further accuracy is required, the method switches automatically to an eigenvalue problem with the augmented matrix. Thus it combines the advantages of the two stages, faster convergence and accuracy, respectively. For the augmented matrix, solving the interior eigenvalue is facilitated by a proper use of the good initial guesses from the first stage and an efficient implementation of the refined projection method. The method can be used with or without preconditioning, on large problems, and can be called with its full functionality from MATLAB through our MEX interface. Numerical experiments illustrate the efficiency and effectiveness of the presented method.