Spacetime geometry with a system of $n$ interacting black holes The Schwarzschild metric is given by $ds^2=-\left(1-\frac{2GM}{rc^2}c^2\right)dt^2+\left(1-\frac{2GM}{rc^2}c^2\right)^{-1}dr^2+r^2d\Omega^2,$ for $d\Omega^2$ the round metric on the $2$-sphere, in local coordinates $(x^0,x^1,x^2,x^3)=(ct,r,\theta,\phi),$ with corresponding metric tensor $$g_{\mu\nu}=\begin{bmatrix} \left(1-\frac{2GM}{rc^2}\right) & 0 & 0& 0\\0&\left(1-\frac{2GM}{rc^2}\right)^{-1}&0&0\\ 0&0&r^2&0\\0&0&0&r^2sin^2\theta\end{bmatrix}.$$

Given this what, then, is the metric (tensor) of a spacetime $(M,\mathcal{O},\mathcal{A},g)$ with a system of $n$ interacting black hole bodies? 

 A: Even Newtonian gravity is extremely complicated if the number of moving bodies exceeds 2. If you add nonlinearity of general relativity and existence of gravitational radiation, then generally, a system of several black holes could be investigated only by numerically solving full system of Einstein equations: a very computationally intensive system of nonlinear PDEs. However, in some cases there are solutions which admit a much simpler analysis. Because of their high degree of symmetry, even in cases where no explicit solutions are known, one can analyze such solutions with much simpler methods. Let us mention some of them. 
Majumdar-Papapetrou  There is a family of Majumdar-Papapetrou (MP) solutions of Einstein-Maxwell equations (general relativity with electromagnetic field). Metric $g$ and  potential $A$ can be
written in the  form
$$\begin{array}{c}
g = -u^{-2}dt^2 + u^2(dx^2+dy^2+dz^2)\,, 
\\
 A = u^{-1}dt\,, 
\end{array}$$
with some nowhere vanishing, say positive, function
$u$ independent of $t$. Einstein-Maxwell equations read then as
$\Delta u=0$, with a 3D flat space Laplacian $\Delta$. A family of solutions for function $u$ has the form
$$
u=1+\sum_{i=1}^I \frac{m_i}{|\vec x - \vec a_i|} \,,
$$
for some positive constants $m_i$. It could be shown that this metric is describing a system of charged black holes with degenerate horizons. Gravitational attraction is exactly compensated by electrostatic repulsion because all charges here are equal to corresponding masses (in appropriate system of units). The original works are:

S.D. Majumdar,
  A class of exact solutions of Einstein’s field equations,
  Phys. Rev. 72 (1947), 390–398.
A. Papapetrou, A static solution of the equations of the gravitational field for an arbitrary charge–distribution, Proc. Roy. Irish Acad.
  A51 (1945), 191–204.

Interpretation as a system of black holes is in following work:

J.B. Hartle, S.W. Hawking,
  Solutions of the Einstein–Maxwell equations
  with many black holes, Commun. Math. Phys.
  26 (1972), 87–101.

But, of course,  there are a lot of more modern treatments, in particular, because such solutions admit Killing spinor and thus could be interpreted as solutions possessing supersymmetry.
Black holes with cosmic strings. There are also static solutions with multiple black holes in which the gravitational attraction is compensated by forces exerted on such black holes by cosmic strings. Cosmic strings exibit $\delta$-like singularity of Ricci tensor and so such solutions could be seen as containing exotic matter. A simplest would be a pair of black holes hanging on from a pair of semi-infinite strings at some distance from each other. Gravitational attraction between black holes is compensated by finetuned string tension.
Black hole universes. Another class of metrics susceptible to simple analysis is "black hole universes" where multiple black holes embedded into cosmological model. For example, one could place black holes in highly symmetric positions on the 3-sphere. Although black holes would be attracting each other, because of the symmetry of their configuration they would be evolving in a cosmological manner with sphere expanding or contracting as a whole. And if we add a finetuned cosmological constant we could even obtain static Einstein universe with multiple black holes! Pioneering work for this would be 

Lindquist, Richard W., and John A. Wheeler. Dynamics of a lattice universe by the Schwarzschild-cell method. Reviews of Modern Physics 29.3 (1957): 432.

And a piece of modern development:

Yoo, Chul-Moon, et al. Black hole universe: construction and analysis of initial data. Physical Review D 86.4 (2012): 044027. arXiv:1204.2411

Higher dimensions. If spacetime dimension $D\ge 5$ then, the no-hair theorem no longer applies. And so among the exact solutions there is such beauty: black saturn. A central spherical black hole is surrounded by a toroidal black ring. Angular momentum of the ring balances its attraction toward the center.

Elvang, H., & Figueras, P. (2007). Black saturn. Journal of High Energy Physics, 2007(05), 050. arXiv:hep-th/0701035.

A: There is no analytic solution of the Einstein equations for a system of more than one black hole. Indeed, the Schwarzschild metric only describes a black hole if:


*

*it is the only object in the universe

*it has already existed for an infinite time and will continue to exist for an infinite time into the future
So given the universe has only existed for 13.8 billion years, and that Hawking radiation means all black holes will eventually evaporate, black holes don't exist.
But leaving this aside, the Schwarzschild metric is a solution to the Einstein equations for a single isolated spherically symmetric mass but it exploits the high symmetry of the situation. If you consider even just two masses no analytic solution can be found. You can't just superpose two Schwarzschild metrics because the Einstein equations are not linear so the sum of their solutions is not a solution.
