We investigate fast and accurate reconstruction methods for susceptibility source separation (S3) in quantitative susceptibility mapping (QSM). S3 separates positive and negative susceptibility sources within a voxel utilizing signal relaxation (R2') for dipole inversion. We propose new primal-dual (PD) methods for S3 and compare them with the alternating Gauss-Newton conjugate gradient (A-GNCG). A-GNCG alters the energy functional, and furthermore its convergence is not guaranteed. In contrast, the proposed PD methods are exact and have convergence guarantees. Validation on a simulated phantom and in-vivo data shows that the PD methods converge faster with better accuracies.