We have implemented a parallel multigrid solver, to solve the initial data problem for 3+1 General Relativity. This involves solution of elliptic equations derived from the Hamiltonian and the momentum constraints. We use the conformal transverse-traceless method of York and collaborators which consists of a conformal decomposition with a scalar $\phi$ that adjusts the metric, and a vector potential $w^i$ that adjusts the longitudinal components of the extrinsic curvature. The constraint equations are then solved for these quantities $\phi$, $w^i$ such that the complete solution fully satisfies the constraints. We apply this technique to compare with theoretical expectations for the spin-orientation- and separation-dependence in the case of spinning interacting (but not orbiting) black holes. We write out a formula for the effect of the spin-spin interaction which includes a result of Wald as well as additional effect due to the rotation of the mass quadrupole moment of a spinning black hole. A subset of these spin-spin effects are confirmed via our numerical calculations, however due to computer time limitations the full parameter space has not yet been surveyed and confirmed. In particular, at the relatively small s eparations ($d \leq 18m$) we are able to consider, we are unable to confirm the expected asymptotic fall-off of $d^{-3}$ for these effects.