In this study, we present an implicit algorithm by applying a primal-dual interior point method (PDIP method) for Gurson-Tvergaard-Needleman model (GTN model) to stabilize the stress update. Although the GTN model widely utilized to realize the change of void volume fraction that dominates ductile fracture in metals, the numerical instability occurs due to the shrinkage of yield surface and the acceleration of void growth. In fact, the implementation of GTN model is limited to the explicit solution method in commercial software such as ABAQUS/Explicit and LS-DYNA. Trial elastic stresses are computed under the assumption that the increment of plastic strain is zero in the conventional return mapping algorithm, which leads to the misjudgment of yield condition since the yield surface shrinks by the evolution of void volume fraction. Thus, it is impossible to judge the yield condition before solving the return mapping equations. In addition, the smoothness of solved equations is required to employ the nonlinear solution method such as the Newton method, whereas the evolution of void volume fraction is approximated as bilinear form to represent the acceleration of void growth. Against these backgrounds, we apply the PDIP method for the stress update of GTN model, in which the inequality constraints in the constitutive model are replaced by an equivalent constrained optimization problem to ensure the numerical stability. After numerical accuracy of the proposed implicit algorithm is verified using iso-error map, its capability is demonstrated throughout several numerical examples that cannot be solved by the conventional return mapping algorithm.