Crystal plasticity theory provides a powerful framework for examining plastic behavior across various scales while integrating crystallographic details into the constitutive model. However, classic rate-independent crystal plasticity models encounter persistent challenges related to the identification of active slip systems and the non-uniqueness of solutions of the underlying set of equations. A widely used and established approach is based on an augmented Lagrangian formulation [1]. Often Lagrangian formulations are solved by means of a fixed-point update that helps to mitigate issues associated with ill-posed crystal plasticity. Although, this fixed-point update is numerically very robust, it is not very efficient. In order to enhance its efficiency, different properties of the algorithm are discussed first. Based thereon, a novel algorithm is proposed that merges the advantages of an augmented Lagrangian approach with nonlinear complementarity problems (NCP). This innovative algorithm also overcomes the ill-conditioning of rate-independent crystal plasticity even without stabilizing hardening or viscous effects. At the same time, it is significantly more efficient than the original augmented Lagrangian formulation. The capabilities and limitations of the proposed algorithm are shown by means of several numerical examples including statistical evaluations of various random crystals as well as RVE simulations of polycrystals.