We present the set of bifurcation tracking algorithms which have been developed in the LOCA software library to work with large scale application codes that use fully coupled Newton's method with iterative linear solvers. Turning point (fold), pitchfork, and Hopf bifurcation tracking algorithms based on Newton's method have been implemented, with particular attention to the scalability to large problem sizes on parallel computers and to the ease of implementation with new application codes. The ease of implementation is accomplished by using block elimination algorithms to solve the Newton iterations of the augmented bifurcation tracking systems. The applicability of such algorithms for large applications is in doubt since the main computational kernel of these routines is the iterative linear solve of the same matrix that is being driven singular by the algorithm. To test the robustness and scalability of these algorithms, the LOCA library has been interfaced with the MPSalsa massively parallel finite element reacting flows code. A bifurcation analysis of an 1.6 Million unknown model of 3D Rayleigh–Bénard convection in a 5 × 5 × 1 box is successfully undertaken, showing that the algorithms can indeed scale to problems of this size while producing solutions of reasonable accuracy.