  We present a novel implicit numerical implementation of the
parabolic-hyperbolic formulation of the constraints of general relativity. The
proposed method is unconditionally stable, has the advantage of not requiring
the imposition of any boundary conditions in the strong field regime, and
offers a holistic (all inclusive) approach to the construction of single and
binary black hole initial data. The new implicit solver is extensively tested
against known exact black hole solutions and is used to construct initial data
for several single and binary black hole configurations.
