The place that you're most likely to find holonomic and non-holonomic systems is in the Langrangian and Hamiltonian formulations of mechanics. The definition of a holonomic system is one whose number of parameters that are required to specify the configuration of the system is equal to the number of degrees of freedom in the system. For a non-holonomic system, there are more parameters than degrees of freedom, where the excess number of parameters is equal to the number of non-integrable equations of constraint. Another way to describe a holonomic system is one whose equations of constraint integrate to functions of just the coordinates, where the constraints are usually functions of the coordinates, time derivatives of the coordinates, etc...The latter definition provides more of a rigorous definition if the intuitive view didn't make sense.
After googling around a bit, it seems that GR is typically viewed as non-holonomic. Analyses of relativistic motion at the quantum level ( ie utilizing quantum potentials and representations ), gets very messy, and starts becoming incredibly difficult to interpret. The closest to an answer I've seen is no-ish ie there are Relativistic systems where the constraints do result in integrable equations, but usually they don't.
As for Phase space. I'll use the Lagrangian example after comparing it to a familiar space. Take a look at $ R^{n} $ for $n=2$. This is the two dimensional space of real numbers. We can geometrically represent this with a graph consisting of two axes. We can label each axis and call it a coordinate. If $ R^{n} $ with it's coordinates, let's use $ x $ and $ y $ , satisfy certain conditions, then we can say that $x$ and $y$ span the space, $R^{n}$. The condition that we have for spanning the space is that $x$ and $y$ axes are linearly independent ie there is no component of the $x$ axis against the $y$ axis and vice versa. More rigorosly: if $ \sum_{i}^{n} \alpha_nx_n = 0 $ for the scalars $ \alpha_n$ not all being 0 where $ x_1 = x $ and $ x_2 = y $. These coordinates generate the space. This is the purely mathematical aspect of phase space.
All a Phase space is, is the space that is spanned by the parameters that specify your system. In the case of Lagrangian mechanics, if I have $ \mathfrak{L}(\dot{q},p) $ then the phase space is the space spanned by taking a the velocity $ \dot{q} $ to be one coordinate for the phase space, and the canonical conjugate momentum, $ p $ to be the other. In this case, the system is specified aptly by these two parameters, so they make up the phase space. These parameters may have explicit dependence on other parameters such as time.
Instead of sticking to standard cartesian coordinates, we generalize our results by writing equations using $q$ as a coordinate. This way, $ q_i $ for $i=1,2,3$, can be cartesian coordinates $x,y,z$ or spherical, cylindrical, curvilinear, or any other kind of coordinates. In this particular case, $ q $ can be regarded as the spatial dimension $x$.
For example, take the harmonic oscillator. What does the phase space for that look like? Let's see!
$$\mathfrak{L} = T - U $$ The kinetic energy, $T$, we take to be $ \frac{1}{2}m\dot{q}^2 $, where the coordinate $ q $ just represents position on a one dimensional line. The potential, $ U $ we take to be $ \frac{kq^2}{2} $. $$ \mathfrak{L} = \frac{1}{2}m\dot{q}^2 - \frac{kq^2}{2} $$ That goes right on into the Euler-Lagrange Equation: $$ \frac{d}{dt}\frac{\partial \mathfrak{L}}{\partial \dot{q}} - \frac{\partial \mathfrak{L}}{\partial q} = 0 $$ Now the holonomic or non-holonomic constraint here that we apply is $ q = \{q:q(t)\} $. This of course means that any time derivatives of $q$ are also functions of some parameter, namely $t$. We take $t$ to be time. Now whether or not the only constraint, and thus the system is holonomic or not is whether or not the differential equation generated by the Euler-Lagrange equation is able to be integrated to a time dependent function, $q(t)$. If we had more dimensions or spatial coordinates, we would have more differential equations, one for movement in each direction as specified by the constraints.
In this case, we get $$ \ddot{q} + \frac{k}{m}q = 0 $$ Differential equations of the form $$ \ddot{q} + \alpha^{2}q = 0 $$ have solutions $$ q = A\cos{\alpha t} + B\sin{\alpha t} $$ Where $A,B=cons.$ are arbitrary constants that result from the initial conditions. This equation holds for any massive particle under the given potential, so let's say solving our initial value problems delivered $ A=1 , B=0 $. Then one exact solution is $$ q = \cos{\omega t}\,\,,\,\,\omega=\sqrt{\frac{k}{m}} $$ So since we were able to determine $q$ as an explicit function of $t$, we have shown that our constraint(s) and thus our system, is holonomic. Since this equation of motion embodies all of the dynamics of the system, we can take a look into the phase space. What are the coordinates of this particular phase space? Whatever the Lagrangian is an explicit function of. These are the parameters that describe our system, since the Lagrangian is a functional of these parameters meant to define said system.
Our Lagrangian was a function of $q$ and $\dot{q}$. Respectively the position and the velocity. I believe for historic reasons, however, we attach the mass term, $m$, to $\dot{q}$ so that our Phase Space is a function of $ p = m\dot{q} $ and $ q $.
Our Phase Space coordinates are functions of a single parameter, so this naturally leads us to being able to create a parametric plot of our phase space of the form $$ < p(t),q(t)>$$ ( Parametric Plots consisting of n functions of 1 parameter generate a curve in n dimensional space )
By calculating the equation of motion for $q$, we just need to solve $$ p = m\dot{q} = -m\omega\sin{\omega t} $$ Giving $$ <-m\omega\sin{\omega t} , \cos{\omega t}>$$ Which plots to
for $ m=\omega=1 $. This should look and feel a bit intuitive an oscillating system creates a circle in phase space. Now since time is our parameter, varying $t$ merely varies where a point is on the circle. If you notice, there is a quantity that stays the same for all $t$: the radius of the circle. More specifically the square of the radius of the circle which turns out to be proportional to the energy of the system! Different radii, different total energies. That of course segues into the concept of the Hamiltonian.
One Lagrangian you might be interested in is $$ \mathfrak{L} = -mc^{2}\sqrt{1 - \frac{\dot{q}^2}{c^2}} $$ This corresponds to a free, relativistic test particle. Pick and chose your constraints ( ie functions to plug in for $\dot{q}$ or its (anti)derivatives. ) And you may be able to find the most specific answer you're looking for. Best of luck, fellow scholar!
