What is the exact gravitational force between two masses including relativistic effects?

Worse than electrodynamics, general relativity is non-linear, in the sense that the field from multiple sources is not just the sum of fields from each isolated source. Even the simple case you are asking about, which is the two-body problem in general relativity, has not been solved exactly.

An even simpler case is the limit $m_2 \to 0$. In that case only $m_1$ affects the geometry of spacetime, and $m_2$ follows a geodesic in that spacetime. This has been solved exactly. The linked article gives details.

[Addendum] To directly answer the original question for the limit $m_2 \to 0$:

Of course if $m_2=0$, then the force between the two masses is $0$. But the thing about classical gravity is that the acceleration due to gravity is independent of mass. (This is the equivalence principle and is actually one of the starting points of the theory of general relativity.) So it still makes sense to ask what would be the acceleration of a negligible mass (aka, a test body) due to the gravity of another mass. The classical answer is $a_2 = \frac{Gm_1}{r^2}$.

In general relativity, the acceleration of a test body due to the gravity of a single spherical, homogenous, non-rotating mass is given exactly by the Schwarzschild solution, which link you can consult for details. The result of which is that

$$\ddot{r} = -\frac{Gm_1}{r^2} + r\dot{\theta}^2 - \frac{3Gm_1}{c^2}\dot{\theta}^2$$ $$\ddot{\theta} = -\frac{2}{r}\dot{r}\dot{\theta}$$ [CORRECTED AND SIMPLIFIED Jan 2]

where $r$ and $\theta$ are polar co-ordinates centred on the gravitating mass, the dots represent differentiation by the proper time of the test body.

So the first term is just the classical radial acceleration $-\frac{Gm_1}{r^2}$. The terms $r\dot{\theta}^2$ and $-\frac{2}{r}\dot{r}\dot{\theta}$ are the classical centrifugal and Coriolis acceleration for polar co-ordinates.

What's not classical is the extra term $\frac{3Gm_1}{c^2}\dot{\theta}^2$. Finally, there is the fact that differentiation is with respect to proper time of the test body. Different test bodies will experience time differently. They can be related by:

$$(1 - \frac{r_s}{r})\dot{t}^2 - \frac{\dot{r}^2}{(1 - \frac{r_s}{r}) c^2} - \frac{r^2\dot{\theta}^2}{c^2} = 1$$

The constant $r_s = \frac{2Gm_1}{c^2}$ is introduced for simplicity. It is called the Schwarzschild radius of the gravitating mass.

Here the co-ordinate $t$ is introduced as a reference time, so $\dot{t}$ is the rate of change of reference time with respect to proper time of the test body. For a distant ($r \to \infty$), stationary ($\dot{r}=0, \dot{\theta}=0$) test body, this becomes $\dot{t} = 1$, so the reference time can be interpreted as the time measured on a distant, stationary clock.

In the classical case of course every body experiences the same time, but it can also be compared to the special relativistic case, where the equation would be:

$$\dot{t}^2 - \frac{\dot{r}^2}{c^2} - \frac{r^2\dot{\theta}^2}{c^2} = 1$$

So what's new in general relativity is the factor $(1 - \frac{r_s}{r})$. To get some idea of scale, for the Earth, $\frac{r_s}{r}$ is about one and a half parts per billion on the surface of the Earth. (Note the Schwarzschild solution is only applicable outside the gravitating body.)

This is exact only for $m_2 \to 0$, but it remains a very good approximation as long as $m_2$ is much smaller than $m_1$, such as for a planet orbiting a star.

Yes the non-relativistic equation of motion

$$\frac{dp}{dt} = F = -\frac{G m_1 m_2}{r^2}$$

for a single particle interacting gravitationally with other is only valid for velocities much smaller than the speed of light and for small enough masses. Notice the minus sign in the expression for the non-relativistic force F --gravitation is attractive--.

First, general relativity is a (geo)metric theory. There is no gravitational forces in general relativity. In general relativity, bodies affected only by gravitation are moving freely, but in a curved spacetime

_{(source: twimg.com)}

The general relativistic equation of motion for a single body is the geodesic equation

$$\frac{DP^\mu}{D\tau} = 0$$

Note the zero at the right, which is a consequence of the absence of gravitational forces in general relativity. The greek indices run over 0,1,2,3 --the spacetime coordinates-- and summation convention is being used. $P^\mu$ is the four-momentum, $\tau$ the proper time and $D$ denotes the covariant derivative, which includes the effects due to spacetime curvature

$$\frac{DP^\mu}{D\tau} = \frac{dP^\mu}{d\tau} + \Gamma_{\nu\lambda}^\mu U^\nu P^\lambda ,$$

where $\Gamma_{\nu\lambda}^\mu$ are the Christoffel symbols and $U^\nu$ the four-velocity.

Second, the field theory of gravity provides a non-geometrical description of gravity. There are gravitational forces in the field theory. In this theory, bodies affected only by gravitation are moving in a flat spacetime --sometimes this is named the flat spacetime approach to gravitation-- but feel a gravitational force associated to gravitons

_{(source: twimg.com)}

The field-theoretic equation of motion for a single body in a gravitational field is the Kalman equation

$$A_\mu^\nu \frac{dP^\mu}{d\tau} = F^\nu = - B_{\mu\lambda}^\nu U^\mu P^\lambda$$

where $F^\nu$ is the exact gravitational force with

$$A_\mu^\nu = \left( 1 - \frac{1}{c^2} \psi_{\lambda\gamma} U^\lambda U^\gamma \right) \eta_\mu^\nu - \frac{2}{c^2} \psi_{\mu\gamma} U^\gamma U^\nu + \frac{2}{c^2} \psi_\mu^\nu$$

and

$$B_{\mu\lambda}^\nu = \frac{2}{c^2} \psi_{\mu,\lambda}^\nu - \frac{1}{c^2} \psi_{\mu\lambda}^{,\nu} - \frac{1}{c^2} \psi_{\mu\lambda,\gamma} U^\gamma U^\nu$$

Here $\psi_{\alpha\beta}$ is the gravitational field potential and the comma denotes the ordinary flat spacetime partial derivative.

As in Electrodynamics, the forces become retarded, the corresponding equations become complicated and include radiation too.