The statement
It seems that the classical programme of the PDE community, i.e., (i) existence (ii) uniqueness (iii) regularity, heavily employing concepts from functional analysis, has not found prominent application in general relativity.
is just plain wrong. You are overlooking quite a lot of stuff. For a modern presentation of the mathematics I would second Ben Whale in recommending Hans Ringstrom's book The Cauchy Problem in General Relativity. Here let me give a bit of historical remarks.
The local wellposedness (after imposing gauge conditions, for reasons already described in Rafe Mazzeo and Peter Michor's answers) of Einstein's equations is a result more than 60 years old! The original proof was published by Yvonne Choquet-Bruhat in 1952 (before she took her current last name), and yes, more modern presentations usually use heavily the notion of Sobolev spaces on manifolds.
Insofar as the local Cauchy problem is concerned, the Einstein system is actually, in some aspects, easier than some of the elliptic and parabolic problems on manifolds. This is principally due to the fact that in appropriate gauge conditions the equations of motion are manifestly hyperbolic (in fact a quasilinear wave equation). Such equations automatically enjoy finite speed of propagation, and thus we can more easily localise the analysis onto individual coordinate patches. That is to say, for the local problem we don't need much of the machinery of global analysis where the geometry and topology of the underlying manifold in the large may come into play.
Perhaps the only conceptually tricky bit of the Cauchy problem for the Einstein system is that, from the very get-go, there is no preferred notion of time on the solution manifold. So whereas in classical non-geometric PDEs the notion of a local existence theorem is stated in the form "there exists $t> 0$ such that a solution exists on $(0,t) \times D$ where $D$ is some domain", the corresponding theorem in general relativity would more naturally look like "there exists a manifold $M$ into which the initial data embeds as a codimension 1 Riemannian manifold with the initial conditions satisfied." This lack of a preferred notion of time also makes it more difficult to interpret what is meant by "global in time solution" to the Cauchy problem.
In 1969 Choquet-Bruhat and Robert Geroch showed that, essentially due to the hyperbolic nature of the equations, one can define the maximal manifold $M$ into which all solutions embed. With this notion one can then formulate the question of "global Cauchy problem" as a question of studying the geometric properties of this maximal solution.
It is well known that this maximal solution can be, in general, quite bad. There are plenty of explicit solutions to Einstein's equation which illustrate this, none-the-least the classical families of black hole space-times. In particular, the classical Schwarzschild solutions are geodesically incomplete in its maximum extension, while the classical Minkowski space is geodesically complete. A natural question, for example, is to ask about the genericity of "geodesic completeness" as a property for maximal solutions of the Cauchy problem, in terms of the initial data. That geodesic incompleteness is generic (in the sense that there exists open sets of initial data that are close to Schwarzschild initial data that leads always to incomplete solutions) turns out to be something that can be established geometrically with minimum amounts of PDE theory. One can read the result off of the incompleteness theorem of Penrose, which establishes a sufficient condition for geodesic incompleteness (the existence of trapped surfaces) and that the existence of trapped surface is an "open condition" and hence is preserved for small perturbations by classical Cauchy stability of the initial value problem at finite times. (A more detailed account is given in Mihalis Dafermos' exposé "The formation of black holes in general relativity (after D. Christodoulou)" for the Seminaire Bourbaki; see Astérisque (2012) 64, Exp. No. 1051.)
The genericity of the case of completeness turns out to be much harder, in terms of the PDE portion. Whereas for the incompleteness case the geometric criterion of Penrose allows us to use "finite time Cauchy stability", which is automatically true in view of the local wellposedness theorems, for the completeness case we have no such criterion and in fact needs to prove global Cauchy stability. In terms of the classical theory of nonlinear PDEs, this corresponds to, roughly speaking, proving global in time existence, with suitable decays, for small data to an initial value problem. Even for simple nonlinear wave equations this is not understood until at least the 70s and 80s. The developments in those two decades led to the understanding that in the physically relevant case of 3 space and 1 time dimensions, the appropriate small data global existence for nonlinear wave equations is not true in general. For Einstein's field equations (which, as described above, can be cast as a quasilinear wave equation, at least locally), this means one needs to search for geometric structure which affords one extra cancellations. As alluded to by Peter Michor, this was done by Demetrius Christodoulou and Sergiu Klainerman in their 500-page opus. (An alternative proof was more recently obtained by Igor Rodnianski and Hans Lindblad; the comparative shorter length testifies to how much the technology for nonlinear wave equations have evolved in the past two decades.)
Along this line of investigation, a further development was produced by Christodoulou in the past few years. As discussed above by Penrose's theorem one automatically has that initial data close to that of Schwarzschild space-time will lead to geodesically incomplete space-times. One may ask whether one can construct data initially far from that of Schwarzschild and still get singularity formation. And the answer turns out to be yes.
The result of Christodoulou and Klainerman on the global stability of Minkowski space actually yielded much more information then just geodesic completeness. It gives precise asymptotic convergence of such a "small data" solution to Minkowski space, and gave a complete description of the "conformal infinity" for such solutions. The aforementioned "quick" proof of stability of geodesic incompleteness by way of Penrose's theorem does not yield nearly as much information. And an ongoing active program of research (involving too many people to be reasonably listed here) is to demonstrate that similar asymptotic stability result with convergence and so forth can be had also for the classical black hole space-times. This is the motivation for the recent bloom in the study of the linear wave equation on black hole space-times.
There are also many other aspects of the global existence problem in general relativity. For example, one can try to consider situations where the initial manifold has compact topology or situations where the initial manifold has symmetries. The former is carefully studied in cosmological settings (I refer again to Ringstrom's book); though oftentimes with high degrees of symmetry so that the system reduces to that of a nonlinear system of ODEs. A spectacular example of the latter is the analysis by Choquet-Bruhat and Vincent Moncrief of solutions with an $U(1)$ symmetry. Under this particular symmetry assumption, and with the so-called CMC gauge condition, the evolution equations of the Einstein vacuum problem reduces to a coupled system of (a) an ODE on the cotangent bundle of Teichmuller space (describing the geometry of the "constant time" slices) and (b) a wave(map) equation.
Coming back to the local well-posedness problem in 3+1 dimensions: in the regime of "classical" solutions, existence, uniqueness, and regularity follows from a very similar analysis to that of Hughes, Kato, Marsden. To close the argument one needs to work in the Sobolev space $H^{5/2+}$ (or classically $H^3$ if you don't want to deal with fractional numbers of derivatives). One may ask about the optimal regularity for the local well-posedness statements. The scaling regularity should be $H^{3/2}$. However, the equation is quasilinear, and it is generally the case that compared to semilinear problems we cannot expect to go all the way down to scaling critical (see for example these two papers of Hans Lindblad 1 2). This problem motivated a body of literature studying the below-classical-regularity local existence problem for quasilinear wave equations. Between the groups of Klainerman-Rodnianski and that of Hart Smith and Daniel Tataru (esp. this paper and this other one), for general quasilinear wave equations the sharp exponent $H^{2+}$ was obtained. For the Einstein vacuum equations specifically, though, one may do better: a recently posted series of pre-prints by Klainerman, Rodnianski, and Jeremie Szeftel obtains a local existence theorem at the regularity level $H^2$ (these are arXiv items: http://arxiv.org/abs/1204.1772 http://arxiv.org/abs/1204.1767 http://arxiv.org/abs/1204.1768 http://arxiv.org/abs/1204.1769 http://arxiv.org/abs/1204.1770 http://arxiv.org/abs/1204.1771 http://arxiv.org/abs/1301.0112).
Some further reading:
